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IMPROVING  THE  NUMERICS  OF  A  THIRD-GENERATION 
WAVE  ACTION  MODEL 


1.  INTRODUCTION 

“SWAN”  (Simulated  WAves  Nearshore)  was  ereated  at  the  Delft  University  of  Teehnology  as  an  exten¬ 
sion  of  third-generation  hyperbolie  wave  aetion  models  to  shallow  water  regions  (e.g.,  Booij  et  al.  1996). 
The  model  has  reeently  beeome  available  in  the  publie  domain,  although  further  model  development  and 
improvement  eontinues.  The  numerieal  seheme  for  spatial  wave  aetion  propagation  used  by  the  model  is  one 
of  several  aspeets  eurrently  being  improved.  Originally,  an  implieit,  first  order,  upwind  numerieal  seheme 
was  ehosen  for  the  model  due  to  the  stability  and  low  eomputational  eost  of  the  seheme.  The  uneonditional 
stability  of  the  seheme  makes  the  model  suitable  for  use  in  shallow  water,  where  eonditionally  stable  models 
(e.g.,  WAM,  WAMDI  Group  1988)  may  beeome  prohibitively  expensive.  The  primary  drawbaek  of  this 
seheme  is  its  relatively  high  rate  of  numerieal  diffusion.  This  diffusion  is  a  major  obstaele  to  the  use  of 
SWAN  for  larger  seale  applieations.  We  have  ineorporated  a  feature  into  an  experimental  version  of  the 
model  that  gives  the  user  the  option  of  using  a  higher  order  and  less  diffusive  numerieal  seheme — a  eyelie 
seheme  taken  from  the  work  of  Stelling  and  Leendertse  (1992). 

2.  MODEL  DESCRIPTION 

Thorough  deseription  and  validation  of  the  wave  model  prior  to  the  addition  of  the  higher  order  numeri¬ 
eal  seheme  ean  be  found  in  the  literature  (Ris  1997,  Booij  et  al.  1999).  Herein,  we  will  only  give  detailed 
attention  to  those  aspeets  of  the  model  that  are  relevant  to  the  geographie  propagation  of  wave  aetion. 

“SWAN”  refers  to  the  offieial  version  v30.75  of  the  eode.  “SWAN-X33”  is  an  unoffieial  extension  of 
v30.75,  modified  by  the  authors.  This  eode  allows  the  use  of  either  the  original  propagation  seheme  or  the 
higher  order  propagation  seheme.  Unless  otherwise  noted,  “SWAN-X”  or  “SWAN-X33”  herein  refers  to  the 
unoffieial  model  with  the  higher  order  seheme  aetivated.  The  ehanges  made  by  the  authors  are  presented  in 
Appendix  A. 

2.1  Governing  Equation:  Action  Balance 

SWAN  is  governed  by  a  two-dimensional  wave  aetion  density  speetrum 


N{x,y,G,6)  -  £'(x,y,(7,0)  /  cr. 


(1) 
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where  x  and  j  denote  geographie  loeation,  a  is  the  relative  frequeney,  and  0  is  the  direetion  of  propagation. 
Wave  aetion  is  propagated  in  geographie  and  speetral  spaee,  while  souree  and  sink  terms  aet  on  the  waves. 
The  aetion  balanee  equation,  in  horizontal  Cartesian  eoordinates  (x,y),  ean  be  written  as 


d  d  d  d  d  S 


(2) 


(e.g.,  Whitham  1974; Phillips  1977;Mei  1983;Hasselmannetal.  1973).  Here, /denotes time, and 5' denotes 
the  total  of  source  and  sink  terms.  The  first  term  represents  the  local  rate  of  change;  the  second  and  third 
terms  represent  geographic  propagation;  the  fourth  term  represents  changes  to  relative  frequency  (e.g.,  by 
nonstationary  depth  or  by  currents);  the  fifth  term  represents  refraction  (by  depth  and  currents).  The  four 
propagation  speeds  are  shown  below  (as  derived  from  linear  wave  theory,  e.g.,  Whitham  1974;  Mei  1983). 

Propagation  speed  in  x-space: 


_dx  _  I 


1  + 


2kd 


sinh2kJ 


ok^ 


2  +V 


(3) 


Propagation  speed  in y- spaee: 


C  =^  =  1 

"  dt  2 


1  + 


2kd 


sinh2kr/ 


oh. 


A 


(4) 


Propagation  speed  in  sigma-space: 


Q  _do  _do 
dt  dd 


—  +  U*Vd 

dt 


j  du 
■  c  k  ^ - 

^  ds  ■ 


(5) 


Propagation  speed  in  theta-spaee: 


"  dt  k 


do  dd  -  du 
——  +  k»  — 
dd  dm  dm 


(6) 


Here,  k  —  {k^,ky)  is  the  wave  number  with  magnitude  k  (related  to  a  through  the  dispersion  relationship  of 

linear  wave  theory),  d  is  water  depth,  U  —  (U^,U  )  is  the  current  velocity,  s  is  the  space  coordinate  in 

direction  0,  and  w  is  a  coordinate  normal  to  s.  The  operator  d/dt  denotes  the  total  derivative  along  a  spatial 
path  of  action  propagation,  and  it  is  defined  as: 


A. 

dt 


(7) 
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where  C^=Yk'  is  the  group  veloeity.  The  term  S  represents  all  souree/sink  terms  that  aet  on  the  aetion 

balanee  equation.  These  terms  inelude  wind  input,  dissipation  (by  whiteeapping,  bottom  frietion,  and  depth- 
indueed  breaking),  and  nonlinear  wave-wave  interaetions  (triads  and  quadruplets).  Detailed  deseription  of 
these  terms  and  how  they  are  implemented  in  SWAN  ean  be  found  in  Ris  (1997)  and  Booij  et  al.  (1999). 

2.2  Treatment  of  Wave  Phase 

The  wave  aetion  balanee  equation  is  a  phase-averaged  formulation.  This  has  the  advantage  of  not  hin¬ 
dering  the  model  with  the  severe  (subwavelength)  resolution  requirements  that  are  typieal  of  phase-resolv¬ 
ing  models.  Unfortunately,  use  of  the  phase-averaged  formulation  means  that  wave  diffraetion  eannot  be 
treated  implieitly  in  the  model  governing  equation.  There  are  tentative  plans  to  add  diffraetion  to  a  later 
version  of  SWAN  in  a  parameterized,  explieit  manner  (this  is  diseussed  further  in  Appendix  B). 

2.3  Numeric  Formulation  of  Propagation 

An  implieit  upwind  differenee  seheme  is  used  for  geographie  and  speetral  propagation  (the  latter  is 
supplemented  with  a  eentered  seheme).  This  seheme,  also  know  as  the  Baekward  Spaee,  Baekward  Time 
(BSBT)  seheme,  has  the  advantage  of  being  uneonditionally  stable  (i.e.,  no  maximum  time  step  size  re¬ 
quired  for  stability).  This  is  advantageous  when  modeling  eoastal  regions  (where  spatial  resolution  is  typi- 
eally  fine),  as  it  ean  be  mueh  more  eeonomieal  than  a  eonditionally  stable  seheme. 

The  finite  differenee  formulation  of  the  governing  equation  is  (from  Booij  et  al.  1999): 


l 

l 

At 

ix  ^ie 

Ax 

Ay 

■(1  +  -  2v[c,nI^  -  (1  -  ■ 

n 

■(1  +  77)hA^].  ,,  -  2ri[c,N].  -  (1  -  il)[ceN].., ' 

2Act 

ix’iy’ie 

2Ae 

where  nisa  time-level  index;  i^,  and  Iq  are  grid  eounters;  and  Ax,  Ay,  Aa,  and  A6  are  inerements  in  time, 
geographie  spaee,  and  speetral  spaee,  respeetively. 


In  speetral  spaee,  the  user  has  the  option  of  supplementing  the  implieit  upwind  seheme  with  an  implieit 
eentered  differenee  seheme.  The  result  is  a  less  diffusive  propagation  in  speetral  spaee,  but  a  greater  likeli¬ 
hood  of  spurious  oseillations.  The  eoeffieients  V  and  r|  are  weighting  faetors,  eontrolling  the  relative  weight 
of  the  two  finite  differenee  sehemes  in  frequeney  and  direetional  spaee,  respeetively. 


2.4  Model  Execution 


2. 4. 1  Nonstationary  Mode 

Running  SWAN  in  nonstationary  mode  is  a  simple  matter  of  time  stepping.  The  user  speeifies  time  step 
inerement,  simulation  begin-time,  and  end-time.  Eaeh  time  step  involves  simple  sweeping  in  geographie 
spaee  and  solution  of  a  matrix  for  A(a,0)  at  eaeh  geographie  grid  loeation.  The  simple  sweeping  in  geo¬ 
graphie  spaee  ean  be  used  with  the  upwind  implieit  numerieal  seheme  beeause  N{x.,y.d^  =  f[N{x.  ^,y.d^, 
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N(x.,y._^,t^,  N(x.,y.,t^_^\  all  three  of  whieh  are  known.  Direetional  spaee  is  divided  into  four  quadrants  that 
are  solved  independently  (exeept  for  interaetions  aeross  the  quadrant  boundaries,  whieh  must  be  aeeounted 
for).  If  eurrents  are  present  or  if  depth  is  not  stationary,  solution  of  a  banded  matrix  is  required.  Otherwise, 
a  simple  tridiagonal  matrix  is  solved.  The  user  has  the  option  of  using  a  set  number  of  iterations  at  eaeh  time 
step;  using  these  iterations,  the  model  more  aeeurately  aeeounts  for  interaetions  aeross  the  internal  bound¬ 
aries  in  speetral  spaee.  Multiple  iterations  at  eaeh  time  step  are  generally  not  used  in  nonstationary  mode, 
under  the  assumption  that  large  ehanges  do  not  oeeur  from  one  time  step  to  the  next. 

The  sequenee  of  operations  in  nonstationary  mode  ean  be  summarized  in  pseudo-eode  as: 

DO  TIME  LOOP 

DO  ITERATION  LOOP  (generally  one  iteration  only) 

DO  QUADRANT  LOOP  (4  times,  onee  per  quadrant/sweep  direetion) 

DO  GEOGRAPHIC  SPACE  LOOPS 

Caleulate  propagation  and  souree  terms. 

Solve  matrix  for  A(a,0). 

(END  LOOPS) 

2.4.2  Stationary  Mode 

The  stationary  mode  of  SWAN  differs  from  the  nonstationary  mode  in  two  respeets.  Firstly,  there  is  no 
time  stepping  (the  solution  is  steady  state).  Seeondly,  iterations  are  required  in  stationary  mode. 

In  stationary  mode,  the  model  is  solved  by  simply  sweeping  through  the  geographie  grid,  with  speetral 
propagation  being  solved  independently  at  eaeh  grid  point.  As  in  nonstationary  mode,  direetional  spaee  is 
divided  into  four  quadrants.  Multiple  iterations  are  neeessary  to  fiilly  propagate  energy  in  speetral  spaee. 
Iterations  allow  transfer  of  energy  aeross  seetor  boundaries  and  faeilitate  the  migration  of  energy  in  fre- 
queney  spaee  due  to  nonlinear  interaetions. 

Convergenee  eriteria  are  used  to  determine  when  to  stop  the  iterations.  Convergenee  eriteria  are  based 
on  a)  relative  and  absolute  ehanges  in  wave  height  between  iterations,  and  b)  relative  and  absolute  ehanges 
in  mean  wave  period.  If  eonvergenee  is  not  reaehed  in  a  user-speeified  maximum  number  of  iterations,  the 
program  terminates. 

2.5  Model  Limitations 

2.5.1  Deficiencies  of  SWAN  at  Oceanic  Scales 

The  SWAN  model  has  a  few  noteworthy  shorteomings,  whieh  tend  to  prohibit  its  use  as  a  general- 
purpose  wave  model,  suitable  for  both  shelf-seale  applieations  O  (10  km)  and  oeean-seale  applieations  O 
(1000  km). 

2.5.1. 1  Diffusion 

The  low  eomputational  eost  and  the  uneonditional  stability  of  the  upwind  implieit  seheme  that  SWAN 
uses  for  geographie  propagation  eomes  at  a  priee — namely,  a  high  degree  of  numerieal  diffusion.  The  term 
“diffusion”  is  used  herein  to  deseribe  the  smoothing  of  wave  aetion  in  spaee.  “Numerieal  diffusion”  refers  to 
diffusion  eaused  by  inaeeuraey  of  a  numerieal  seheme.  “Physieal  diffusion”  does  not  exist  in  SWAN,  as  the 
governing  equations  are  equations  of  adveetion  rather  than  adveetion-diffusion;  thus  all  diffusion  in  SWAN 
is  numerieal.  In  this  report,  we  are  primarily  eoneemed  with  the  numerieal  diffusion  by  SWAN  in  geo- 
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graphic  space  (as  opposed  to  spectral  space).  It  should  be  noted  that  the  numerical  diffusion  in  SWAN  is  a 
mass-conserving  process:  wave  action  is  not  lost  through  diffusion,  except  in  cases  where  wave  action 
moves  out  of  the  system  through  the  open  boundaries  as  a  result  of  diffusion. 

Figure  1  shows  a  severe  example  of  this  numerical  diffusion.  As  evidenced  by  the  scale  of  Fig.  1,  this 
diffusion  is  not,  strictly  speaking,  related  to  scale.  However,  two  factors  make  diffusion  a  greater  concern  at 
oceanic  scales  than  at  shelf  scales: 

1 .  Diffusion  is  exacerbated  by  nonstationary  input  conditions,  which  are  more  typical  for  larger-scale 
applications. 

2.  Grid  point-to-grid  point  variations  in  wave  energy,  either  locally  generated  or  specified  at  the 
boundaries,  are  typically  greater  for  larger  scale  applications.  This  results  in  a  higher  rate  of  diffu¬ 
sion. 


X  (dimensionless) 

Fig.  1  —  Idealized  test  ease:  propagation  of  wave  energy  past  a  submerged  breakwater  with  a  gap 
using  the  diffusive  first  order  upwind  seheme.  Wave  energy  is  shown.  Energy  is  speeified  at  the  left 
(x  =  0)  boundary:  energy  =  2  at  gap,  energy  =  1  elsewhere.  Energy  is  propagated  at  an  angle  with  the 
grid.  Model  parameters:  CEL  =  1.4;  =  8.0  (dimensionless);  dx  =  dy  =  0.5  (dimensionless);  diree- 

tion  of  transport  theta  =  35  degrees  from  v  axis;  and  201  grid  points  in  v;  and  201  grid  points  iny. 
Distanee  is  nondimensionalized  by  Courant  number. 


Increasing  the  applicability  of  SWAN  by  reducing  this  numerical  diffusion  is  the  primary  focus  of  this 
report. 

2.5.1.2  Nonstationary  Input  at  Boundaries 

Nonstationary  input  at  boundaries  has  not  yet  been  implemented  in  SWAN  (except  in  the  context  of 
nested  simulations).  At  the  large  time  scales  of  oceanic  applications,  this  is  a  rather  severe  omission.  Cur¬ 
rently,  there  are  plans  to  implement  nonstationary  boundary  input  in  a  public  version  to  be  released  in  1999. 
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2.5.1.3  Spherical  Coordinates 

The  curvature  of  the  Earth  becomes  important  at  larger  geographic  scales.  Therefore,  ocean-scale  wave 
models  (e.g.,  WAMDI  1988)  are  typically  calculated  in  spherical  coordinates.  The  present  version  of  SWAN 
does  not  allow  for  spherical  coordinates,  although  there  are  plans  for  a  future  version  that  will. 

2.5.2  Other  Limitations 

2.5.2. 1  Diffraction 

As  stated  in  Section  2.2,  diffraction  cannot  be  included  implicitly  in  the  governing  equation  of  a  hyper¬ 
bolic  wave  model  such  as  SWAN.  The  developers  of  SWAN  are,  however,  investigating  the  possibility  of 
including  diffraction  in  parameterized  form  in  a  later  release  of  SWAN. 

2.5.2.2  Open  Boundaries 

Boundary  conditions  in  both  geographic  and  spectral  space  are  fiilly  absorbing.  Wave  energy  enters  the 
domain  along  the  open  boundaries  only  if  specified  by  the  user.  Thus,  when  the  lateral  boundaries  are  not 
specified,  this  generally  results  in  a  shadow-region  adjacent  to  one  or  both  of  the  lateral  boundaries. 

Awareness  of  this  feature  is  especially  important  for  coastal  simulations.  In  a  typical  coastal  simulation, 
the  model  user  prescribes  wave  conditions  along  only  the  offshore  boundary  of  a  rectangular  grid  domain 
(e.g..  Fig.  2).  In  order  to  prevent  the  shadow-region  from  influencing  the  region  of  interest  (which  is  often 
close  to  shore),  the  user  must  increase  the  lateral  extent  of  the  domain.  Specifying  wide  directional  spectra 
(e.g.,  -80°  to  +80°)  as  model  input  along  a  single  open  boundary  would  require  the  use  of  a  domain  with  a 
wide  aspect  ratio. 

2.5.2.3  Hydrodynamics 

SWAN  does  not  attempt  to  quantify  the  effect  of  wave  climate  on  circulation  in  the  domain.  This  is  well 
beyond  the  scope  of  the  model.  Currents  are  present  as  input  only. 

3.  CHOOSING  THE  NEW  SCHEME 

The  purpose  of  this  study  is  to  reduce  the  artificial  numerical  diffusion  in  SWAN  by  replacing  the 
diffusive,  upwind,  implicit  scheme  (used  for  geographic  propagation  of  wave  energy)  with  a  less  diffusive, 
higher-order  numerical  scheme.  To  this  end,  we  have  evaluated  the  capabilities  of  various  potential  replace¬ 
ment  schemes.  First,  we  compared  a  range  of  schemes  in  one-dimensional  mode.  Then  we  conducted  two- 
dimensional  tests  using  several  of  the  more  promising  schemes. 

3.1  Tests  of  One-Dimensional  Propagation 

3.1.1  Description  of  One-Dimensional  Tests 

Two  basic  types  of  one-dimensional  tests  were  performed: 

1.  A  spike  of  wave  energy  defined  as  the  initial  condition  of  the  model: 


N{x)  =  {64  X  [(V  -  0.5)^  -ye,]fif%<x<%, 
N{x)  =  0  otherwise. 


(9) 
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Fig.  2  — An  example  of  the  shadow  region  ereated  if  wave  energy  is  not  speeified  along  the 
lateral  boundaries.  Waveheight,  is  shown.  Here,  input  is  speeified  at  the  left  (x  =  0)  boundary, 
and  is  not  speeified  at  the  (y  =  0)  boundary  or  the  (y  =  40  km)  boundary.  The  (x  =  40  km) 
boundary  lies  on  land.  There  is  a  gradual  transition  between  the  illuminated  region  and  shadow 
region  due  to  the  broadness  of  the  input  direetional  speetra  and,  to  a  lesser  extent,  numerieal 
diffusion. 


A  single  group  velocity,  C  ,  is  specified  and  the  spike  is  then  propagated  by  time  stepping.  N  at  the  two 
boundaries  is  set  to  zero  for  all  time  steps. 

2.  A  sinusoidal  variation  of  wave  energy  is  defined  as  an  unsteady  open  boundary  condition: 


=  1  +  0.1  xsin(co0. 


(10) 


The  initial  condition  is  a  state  of  rest  {N(x)=l  for  all  x). 

3.1.2  Scheme  Descriptions  and  Observations  of  Performance 

The  numerical  schemes  are  described  in  the  following  paragraphs.  The  advantages  and  disadvantages 
of  each  scheme  are  listed.  These  are  based  on  preliminary  knowledge  of  accuracy  and  stability  conditions,  as 
well  as  observations  of  one-dimensional  test  results. 

Where  possible,  the  finite  difference  approximations  of  spatial  and  temporal  derivatives  are  given  sepa¬ 
rately.  In  the  context  of  the  one-dimensional  tests,  these  approximations  are  used  in  the  one-dimensional 
hyperbolic  simple  wave  equation: 
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dt  ^  dx 


0. 


(11) 


Stability  requirements  stated  in  this  seetion  are  only  applieable  to  the  one-dimensional  simple  wave 
equation. 

3. 1.2.1  BSBT  (Backward  Space,  Backward  Time) 

Also  known  as  the  first-order,  upwind  implieit  seheme,  this  is  the  numerieal  seheme  used  for  geo- 
graphie  propagation  in  SWAN.  Approximations  to  the  derivatives  using  this  seheme  are: 


dN 

dt  At 

and 

37V  ^  -  Nl, 

dx  Ax 


(12) 


(13) 


Advantages: 

•  Compaet:  requires  only  three  points  (ideal  for  use  at  boundaries) 

•  Fast  eomputationally 

•  Uneonditionally  stable 

Disadvantages: 

•  Severe  numerieal  diffusion 

3. 1.2.2  Second-Order,  Upwind,  Implicit 

This  scheme  is  the  second-order  extension  of  the  BSBT  scheme  (e.g.,  Roache  1972).  Applied  to  the 
one-dimensional  simple  wave  equation,  the  scheme  gives: 


^  -  1^^(^  +  1)(7V;  -  2Nl,  +  Nl^). 


(14) 


Advantages'. 

•  Unconditionally  stable 

•  Less  diffusive  than  BSBT  scheme 

Disadvantages’. 

•  High  rate  of  diffusion  relative  to  most  other  higher-order  schemes 

•  Unlike  the  first-order  scheme,  creates  nonphysical  oscillations 
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3.1.2.3  SLl:  Stelling  and  Leendertse  (q^=  0,  q  =  1/6) 

This  scheme  is  taken  from  the  work  of  Stelling  and  Leendertse  (1992),  using  the  parameters  qg=0  and 
q^=\l6  in  the  derivation.  The  scheme  uses  a  four-point  upwind  differencing  at  time  level  t^  and  a  three  point 
centered  differencing  at  time  level  t^ 


dN  _ 

dt  At 


(15) 


and 


dN  ION" -ISN-.+eN-.-N"  A 

_  ^  l-rl _ l  —  j  I  _ I _ l  —  j _ t  — Z _ l  —  J 

dx  2 2Ax  6Ax  ^ 


Advantages: 

•  Unconditionally  stable 

•  Less  diffusive  than  BSBT  scheme 


(16) 


Disadvantages: 

•  Some  numerical  diffusion  (increasing  with  higher  Courant  numbers  and  lower  resolution).  At  higher 
Courant  numbers,  this  diffusion  is  large  relative  to  most  other  schemes  tested. 

•  Some  nonphysical  oscillations. 

Figure  3  shows  example  comparisons  of  BSBT  and  SLl  solutions. 

3. 1.2.4  SL2:  Stelling  and  Leendertse  {q  =  0,  q  =  0) 

This  scheme  employs  a  three-point  upwind  differencing  at  time  level  t^  and  a  three-point  centered 
differencing  at  time  level  t^_^: 


and 


dN  _ 

dt  At 


dx 


A/'"”* 


i+2 


+ 4N";'  -  4N":'  -  sn"  -  4a^;_,  + 


12Ax 


2Ax 


Advantages'. 

•  Unconditionally  stable 

•  Fewer  nonphysical  oscillations  than  SLl  at  higher  Courant  numbers 
Disadvantages'. 


(17) 


(18) 


Some  numerical  diffusion  (more  than  the  SLl  scheme) 
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Fig.  3  —  A  one-dimensional  sinusoidal  test  ease  eomparing  the  BSBT  and  SLl  sehemes.  Model  parameters: 
CFL=  0.5;  =  8.0  (dimensionless);  Ar  =  4  (dimensionless);  20  points  per  wavelength. 


3.1.2.5  SL3:  Stelling  and  Leendertse  {q  =  -1/6,  q  =  0) 

This  scheme  employs  a  three-point  upwind  differeneing  at  time  level  and  a  five-point  eentered 
differeneing  at  time  level 


dN  _ 

dt  At 

and 

_  1  ( -  v;_~‘  3v;-4v;_,  +  v;_,^ 

dx  2 2Ax  2Ax  ^ 

Advantages’. 

•  Uneonditionally  stable 

•  Fewer  nonphysieal  oseillations  than  SLl  at  higher  Courant  numbers 
Disadvantages’. 

•  Some  numerieal  diffusion  (more  than  the  SLl  seheme) 


(19) 


(20) 


Figure  4  shows  some  example  eomparisons  of  the  Stelling  and  Leendertse  sehemes. 
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Fig.  4 - A  one-dimensional  spike  propagation  test  ease  eomparing  the  three  Stelling  and  Leendertse  sehemes. 

Model  parameters:  CFL  =  1.6;  =  8  (dimensionless);  Ax  =  1/32  (dimensionless). 


3.1.2.6  Box 


The  Box  scheme  (e.g.,  Petit  1997)  is  a  compact  scheme  which  uses  a  two-point  upwind  differencing  at 
both  time  level  t  and  at  time  level  /  , : 

n  n-\ 


and 


dt  2 


m  -  Nr'  K ,  -  Nr 

i  i  z  — 1  z  — 1 


n-\  \ 


At 


At 


dx  2 


x;  -  Nr  n; 


n-\ 


Ax 


+  ■ 


■N' 


n-l  \ 


i-l 


Ax 


(21) 


(22) 


Advantages: 

•  Unconditionally  stable 

•  Good  diffusion  characteristics  at  all  resolutions;  in  general,  less  diffusive  than  the  SLl  scheme 

•  Compact  (ideal  for  use  at  boundaries) 

Disadvantages: 

•  Relatively  high  amount  of  nonphysical  oscillation,  especially  at  higher  Courant  numbers 

•  Phase  speed  error  is  significant  at  lower  Courant  numbers 
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3.1.2.7  Theta-Box 

Theta-Box  (e.g.,  Petit  1997)  is  a  generalized  form  of  the  Box  seheme,  weighting  the  {t^  differeneing 
and  the  differeneing  in  the  spatial  derivative  using  a  eoeffieient  0: 


and 


dt 


N"  -  Nl,  -  ^ 


St 


■  + 


dN  if  N""  -N""  ^ 

/Vi  yv,_i 


dx  2 


Ax 


Ax 


J 


(23) 


(24) 


Advantages: 

•  In  general,  less  diffusive  than  the  SLl  seheme 

•  Compaet  (ideal  for  use  at  boundaries) 

•  Uneonditionally  stable  for  6  >  0.5 

Disadvantages: 

•  Relatively  high  amount  of  nonphysieal  oseillation 

•  At  0  0.5,  this  seheme  generally  performed  worse  than  the  standard  Box  seheme  (0  =  0.5) 

3.1.2.8  Crank-Nicholson 

The  Crank-Nieholson  seheme  (e.g.,  Roaehe  1972)  is  a  eomputationally  implieit  seheme  {N  values  at 
time  level  ^^eannot  be  solved  independently),  the  eomponents  of  whieh  are: 


dN  _ 

dt  At 

and 

dN  __  1  ( yv;^i  -  yv"_i  ^  n;-'  -  v;_~‘  ^ 
dx  2 2  Av  2Ax  ^ 

Advantages: 

•  Uneonditionally  stable 
Disadvantages: 

•  Very  slow.  Beeause  Crank-Nieholson  is  a  eomputationally  implieit  seheme,  a  matrix  solver  is  re¬ 
quired. 

•  Moderate  amount  of  diffusion 

•  Large  amount  of  nonphysieal  oseillation 


(25) 

(26) 
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3.1.2.9  Borsboom 

The  Borsboom  scheme  (e.g.,  Petit  1997)  is  another  computationally  implicit  scheme.  It  is  similar  to  the 
Crank-Nicholson  scheme,  with  a  time  derivative  that  is  spatially  weighted: 


and 


dt 


1^ 
61 


N’' 


■N: 


i+l 


At 


4 

+  - 

6 


N"  -  NT 


At 


1^ 
+  - 

6 


NU-N[ 


n-\  \ 


V 


At 


dx 


_  M' 


n-l  \ 


2Ax 


2Ax 


(27) 


(28) 


Advantages'. 

•  Unconditionally  stable 

•  Less  diffusive  than  Crank-Nicholson  or  SLl 

Disadvantages'. 

•  Very  slow  (in  terms  of  CPU  time,  it  is  the  most  expensive  of  all  schemes  tested  here) 

•  Moderate  nonphysical  oscillation 

Figure  5  compares  the  SLl  scheme  to  the  two  computationally  implicit  schemes. 

3.1.2.10  NISL:  Non-Interpolating  Semi-Lagrangian 

This  scheme,  from  Olim  ( 1 994),  is  a  modification  of  the  Lax- Wendroff  method.  It  uses  a  finite  differencing 
based  on  three  centered  upwind  points  at  a  distance  in  the  grid  determined  by  the  Courant  number.  At  time 
level  the  scheme  “looks  back”  along  the  characteristic  line  to  three  neighboring  points  at  time  level  t^_^. 
The  application  of  the  scheme  to  the  one-dimensional  simple  wave  equation  is  written  as: 


n;  -  n::1  - - pXnt^  -  + y2(r - pfiK_%  - in::1  - 


(29) 


where  p  is  the  Courant  number,  and  p  is  the  nearest  integer  to  |LL. 

Advantages'. 

•  Unconditionally  stable 

•  No  inherent  upper  limit  to  Courant  number 

•  Exceptional  performance  at  higher  CFLs  (e.g.,  at  |LL  >  2,  performs  better  than  any  other  scheme 
tested) 
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Fig.  5  — A  one-dimensional  spike  propagation  test  ease  eomparing  the  SLl  seheme  to  the  two  eomputationally 
implieit  sehemes.  Model  parameters:  CFL  =  1.5;  =  1  (dimensionless);  Ax  =  1/32  (dimensionless). 


Disadvantages: 

•  At  lower  Courant  numbers,  exhibits  signifieant  diffusion  and  nonphysieal  oseillation  (eompared  to 
eomparable  sehemes) 

•  The  numerie  steneil  (the  points  that  the  seheme  uses)  varies  aeeording  to  Courant  number,  whieh 
eould  be  ineonvenient  when  eoding  a  model 

3.1.2.11  BSFT  (Backward  Space,  Forward  Time) 

This  first-order,  upwind,  explieit  seheme  is  stated  as  (e.g.,  Roaehe  1972): 


and 


dN 

dt  At 

dN  _ 

dx  Ax 


(30) 


(31) 


This  is  the  primary  numerical  scheme  used  in  the  third-generation  hyperbolic  wave  model  WAM  (WAMDI 
Group  1988). 
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Advantages: 

•  Compact:  requires  only  three  points 

•  Exceptionally  fast  computationally  (the  fastest  of  all  the  schemes  tested  here) 
Disadvantages: 

•  Severe  numerical  diffusion  (although  generally  less  than  BSBT  scheme) 

•  Conditionally  stable  (stable  for  Courant  number,  q  <  1.0) 

•  Generally  performs  well  only  when  Courant  number  is  one  or  slightly  less  than  one 

3.1.2.12  Lax 

The  Lax  scheme  is  (Roache  1972): 


and 


dt 


2 


V 


At 


dN  ^ 

dx  2Ax 


(32) 


(33) 


Advantages: 

•  Fast  computationally 
Disadvantages: 

•  Conditionally  stable 

•  High  diffusion,  especially  at  lower  Courant  numbers 

•  Worst  performer  of  all  explicit  schemes  tested:  generally,  performs  well  only  when  Courant  number 
is  one  or  slightly  less  than  one 

3.1.12.13  BCMLn3:  Backward  Characteristic  Method,  Lagrangian  Interpolation  {n  =  3) 

The  formulation  of  this  method  is  similar  to  that  of  NISL  in  that  it  tracks  the  solution  backward  along 
characteristic  rays.  The  BCMLn3  scheme  interpolates  to  the  surrounding  points  at  the  end  of  the  backward 
projection  using  Lagrangian  interpolation  functions.  The  approximation  of  the  one-dimensional  hyperbolic 
wave  equation  using  this  method  is  (e.g..  Petit  1997): 

2)(^  -  + 1)(^ + 2)Ni-,^ 

3)iil  - 

+  Mr-  3)(jU  -  + 1)(^  +  2) (34) 

-Mr  -  -  m + m + 2)ivr‘ 

+Mr-  + 2)A^r+i‘ 

-Mr-  Mr  -  Mr  -  Mr)(r + m:;2 
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Advantages: 

•  At  lower  Courant  numbers  (e.g.,  |LL  <  2),  is  generally  more  aeeurate  than  any  other  seheme  tested 
here 

Disadvantages: 

•  Conditionally  stable:  |LL  <  1  (ealeulated);  |Li  <  2  (observed) 

3.1.2.14  Fromm 

The  Fromm  seheme  (Fromm  1968)  uses  one  point  at  time  level  {t^  and  four  points  at  time  level  {t^  ^). 
The  one-dimensional  wave  equation  is: 


n;  «  A^;-‘  -  KMC/  -  Ci‘)  +  -  2A^r‘  +  Ci‘) 

11(11-  l)(C“i‘  -  3A^r‘  +  3Ci‘  -  K-A 


Advantages: 

•  Generally  less  diffusive  than  SLl  seheme  (although  more  diffusive  than  BCMLnS) 

•  Zero  phase  error 

Disadvantages: 

•  Conditionally  stable:  |LL  <  1 

3.1.2.15  QUICKEST:  Backward  Characteristic  Method,  Lagrangian  Interpolation  {n  =  2) 

This  scheme  is  stated  as  (Leonard  1979): 

N" «  A^;-‘  - 'A /i(2N.~'  + 

+  KM(Ci‘  -  +  NI\')  +  Aii\N.:2  -  3N.:'  + 


This  scheme  is  used  by  the  third-generation  hyperbolic  wave  model  WAVEWATCH  (Tolman  1995). 
Advantages: 

•  Generally  less  diffusive  than  SLl  scheme  (although  more  diffusive  than  BCMLnS) 
Disadvantages: 

•  Conditionally  stable:  p  <  1  (calculated);  p  <  2  (observed) 

Figure  6  shows  an  example  of  a  comparison  of  the  SLl  to  the  five  conditionally  stable  schemes  tested, 
including  QUICKEST  In  this  case,  and  in  most  other  cases,  the  higher-order  conditionally  stable  schemes 
(QUICKEST,  Fromm,  and  BCMLnS)  are  more  accurate  than  SLL 
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Fig.  6  - A  one-dimensional  spike  propagation  test  ease  eomparing  the  SLl  seheme  to  the  five  eonditionally 

stable  sehemes.  Model  parameters:  CFL  =  0.75;  Cg  =  I  (dimensionless);  Ax  =  1/32  (dimensionless). 


3.2  Tests  of  Two-Dimensional  Propagation 

After  the  one-dimensional  tests  were  eondueted,  most  of  the  sehemes  were  eliminated  as  eandidates  for 
the  new  model.  In  order  to  preserve  the  uneonditional  stability  of  the  model,  the  eonditionally  stable  sehemes 
were  eliminated.  Due  to  eonsideration  of  eomputational  eosts,  the  two  eomputationally  implieit  sehemes 
were  eliminated.  The  theta-Box  seheme  was  also  eliminated,  as  it  offered  little  apparent  advantage  over  the 
Box  seheme.  The  seeond-order,  upwind,  implieit  seheme  was  eliminated  due  to  poor  performanee  in  the 
one-dimensional  tests.  The  following  sehemes  were,  therefore,  tested  in  two  dimensions:  BSBT,  SLl,  SL2, 
SL3,  Box,  and  NISL.  The  SL2  and  SL3  sehemes  were  ineluded  in  only  two  of  the  two-dimensional  tests. 
These  two  sehemes  were  subsequently  eliminated,  as  performanee  was  eonsistently  inferior  to  that  of  the 
SLl  seheme. 


The  NISL  seheme,  in  two  dimensions,  is  stated  by  Olim  (1994)  as: 


N""  (1-  I 

h^rx  ^  y^^ i-p-lj-q-l  ^  i-p+lj-q+l  i-p+lj-q-l  i-p-lJ-q+G 


(37) 


where  is  the  residual  Courant  number  inx-spaee,  (Q  —  p);  pis  the  nearest  integer  to  \i^=  (C^ 

is  the  residual  Courant  number  iny-spaee,  (C^  —  ^);and  q  is  the  nearest  integer  to  |Li^=  (C^  "%>')• 


The  Box  seheme  was  extended  to  two-dimensional  spaee  by  rewriting  it  as  a  “eube”  seheme,  where 
eaeh  derivative  is  ealeulated  from  the  eight  points  on  an  upwind-oriented  eube: 
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and 


dt  4 


-  A/^"“^  A/^"  -  A/^”“^ 


V 


At 


■  + 


Ar 


Ar 


+  ■ 


n-\  \ 

j 


At 


(38) 


^~i 

3x  4 


W"  _  W"  W"-l  _  W"-l  _  M"  W"-l  _  W«-l  ^ 


Ax 


Ax 


Ax 


Ax 


(39) 


dN  1 
4 


^  A/-:‘ 


V 


Ay 


Ay 


j^n-l  j^n  _j^n  j^n-\  _  J^n-l 


Ay 


Ay 


(40) 


In  the  cases  of  the  other  four  schemes  (BSBT,  SLl,  SL2,  SL3),  the  equations  for  partial  derivatives 
given  in  Section  3.1.2  were  applied  directly  to  the  two-dimensional  simple  wave  equation: 


dt  "  dx  "  dy 


(41) 


3.2.1  Nonstationary,  Nonuniform  Input 

3.2. 1.1  Two-Dimensional  Propagation  of  a  Spike 

A  spike  of  the  following  form  was  propagated  using  the  six  models  described  above: 

N{x,  y)  =  [cosh(v/'  *  T)]"' ,  (42) 


where  \|/  is  a  narrowness  coefficient,  and  T  is  the  distance  from  the  spike  center.  Examples  of  results  are 
shown  in  Figs.  7  through  10.  Several  observations  were  made: 

•  As  expected,  relative  performance  of  the  schemes  was  a  function  of  Courant  number,  and  (to  a 
lesser  degree)  spatial  resolution.  The  NISL  scheme  is  clearly  superior  at  larger  Courant  numbers 
(e.g.,  p  >  1),  while  the  Box  and  the  Stelling  and  Leendertse  schemes  do  best  at  lower  Courant 
numbers  (e.g.,  p  <  0.5). 

•  As  in  the  one-dimensional  tests,  the  SLl  scheme  generally  creates  less  diffusion  than  the  other  two 
Stelling  and  Leendertse  schemes. 

•  SLl  generally  has  slightly  greater  diffusion  than  the  Box  scheme. 

•  Box  exhibited  more  problems  with  nonphysical  oscillation  than  SLl . 

•  The  NISL  scheme  sometimes  becomes  unstable  in  cases  where  (p^^~  p^^~  0.5).  Explanation  for  this 
phenomenon  could  not  be  found  via  stability  analysis. 
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Fig.  7  —  Results  of  two-dimensional  spike  propagation  test  using  the  BSBT  seheme.  The 
initial  shape  and  loeation  of  the  four  spikes  is  shown  in  the  lower  left  eomer.  The  four  spikes 
are  propagated  at  angles  of  0,  15,  30,  and  45  degrees.  The  parameter  “D”  (shown  above  the 
end  loeation  of  eaeh  of  the  four  spikes)  indieates  the  fraetion  of  spike  amplitude  retained,  a 
rough  measure  of  diffusion  (lower  D  ==  greater  diffusion).  Model  parameters:  CFL  =  0.75;  Ax 
=  =  0.1  (dimensionless);  =  0.18  (dimensionless). 
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Fig.  8  — Results  of  two-dimensional  spike  propagation  test  using  the  Box  seheme. 
Model  setup  is  identieal  to  that  of  the  BSBT  test  ease  deseribed  in  Fig.  7. 
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Fig.  9  —  Results  of  two-dimensional  spike  propagation  test  using  the  NISL  seheme. 
Model  setup  is  identieal  to  that  of  the  BSBT  test  ease  deseribed  in  Fig.  7. 


-0.8 


-0.6 


0.4 


0.2 


-0.2 


5  10  15  20  25 

X  (dimensionless) 

Fig.  10  — Results  of  two-dimensional  spike  propagation  test  using  the  SLl  seheme. 
Model  setup  is  identieal  to  that  of  the  BSBT  test  ease  deseribed  in  Fig.  7. 
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3.2. 1.2  Sine  Wave  Propagation  Test 

Two-dimensional  “waves  of  wave  energy”  were  simulated  in  these  tests.  Conditions  were  speeified 
along  the  {x  =  0)  and  (y  =  0)  boundaries: 

(t)  =  1  -0  +  sin(3;  xk^-axt)  (43) 

and 

(x)  =  1.0  +  sin(x  xk^-GXt).  (44) 


Figure  11  shows  the  results  from  one  of  these  simulations,  and  Figs.  12  and  13  show  transeets  of  wave 
energy  along  the  direetion  of  energy  propagation. 
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Fig.  11  —  Example  result  of  the  two-dimensional  sine  wave  propagation  test  for  the  Box  seheme. 
Global  wave  aetion  is  shown.  Model  parameters:  CFL  =  0.75;  direetion  of  transport  theta  =  30 
degrees  from  x  axis;  15  points  per  wavelength. 


Unlike  the  spike  propagation  tests,  these  tests  showed  a  marked  differenee  in  diffusion  eharaeteristies 
between  the  Box  and  SLl  sehemes.  The  SLl  seheme  often  exhibited  severe  diffusion,  espeeially  at  higher 
Courant  numbers,  whereas  diffusion  levels  were  generally  low  with  the  Box  seheme. 

Just  as  in  the  one-dimensional  tests,  the  SL2  and  SL3  sehemes  exhibited  greater  diffusion  than  the  SLl 
seheme,  while  offering  only  a  minor  deerease  in  nonphysieal  oseillations. 
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Fig.  12  — A  two-dimensional  sine  wave  propagation  test:  wave  aetion  along  the  direetion  of  propagation. 
Model  parameters:  CFL  =  0.25;  direetion  of  transport  theta  =  45  degrees  from  v  axis;  15  points  per  wavelength. 
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Fig.  13  —  A  two-dimensional  sine  wave  propagation  test:  wave  aetion  along  the  direetion  of  propagation. 
Model  parameters:  CFL  =  1.0;  direetion  of  transport  theta  =  30  degrees  from  v  axis;  15  points  per  wavelength. 
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The  sine  wave  tests  provided  greater  insight  into  errors  in  propagation  speed  eaused  by  inadequate 
spatial  resolution,  unfavorable  Courant  numbers,  or  both.  Observations  of  nonphysieal  oseillation  showed  a 
similar  trend  to  that  of  diffusion:  the  NISL  seheme  exeels  at  higher  Courant  numbers,  while  Box  and  SLl  do 
best  at  lower  numbers. 

3.2.2  Stationary,  Uniform  Input:  Test  of  Leading  Oscillations 

This  is  a  very  simple  test  with  eonstant  and  uniform  input  along  the  offshore  (x  =  0)  boundary.  The 
initial  eondition  is  a  state  of  zero  energy  throughout  the  grid.  Thus,  the  test  is  similar  to  the  shoek  propaga¬ 
tion  tests  often  used  for  eomputational  fluid  meehanies  eodes.  In  order  to  study  oseillations  at  the  leading 
edge  of  the  energy  front,  the  simulations  were  not  run  to  steady  state.  Results  were  eompared  along  transeets 
in  the  direetion  of  energy  propagation.  Figures  14  and  15  show  some  of  the  results  of  these  tests.  Beeause  of 
the  emphasis  on  nonphysieal  oseillation,  observations  were  very  favorable  toward  the  SLl  seheme.  The 
NISL  seheme  and  Box  seheme  exhibited  severe  oseillations,  with  no  signifieant  advantage  over  the  SLl 
seheme  at  any  spatial  resolution  or  time  step  inerement. 


Fig.  14 — A  two-dimensional  test  of  leading  oseillations:  wave  aetion  along  the  direetion  of 
propagation.  Model  parameters:  CFL  =  0.25;  direetion  of  transport  theta  =  45  degrees  from 
V  axis;  Av  =  Ay  =  0.1  (dimensionless). 


3.2.3  Stationary,  Nonuniform  Input:  Submerged  Breakwater  Tests 

A  submerged  breakwater  was  erudely  simulated  by  speeifying  wave  heights  at  the  offshore  boundary: 

H^fx  =  0)  =  2.0  m,  where  4700  m  <  y  <  5 100  m,  and 
H^/x  =  0)  =  1.0  m,  where  y  <  4700  m  andy  >  5100  m. 
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Fig.  15 - A  two-dimensional  test  of  leading  oseillations:  wave  aetion  along  the  direetion 

of  propagation.  Model  parameters:  CFL  =  1.2;  direetion  of  transport  theta  =  45  degrees  from 
V  axis;  Ax  =  =  0.1  (dimensionless). 


This  input  was  stationary  (constant  in  time),  but  the  simulations  were  not  run  to  steady  state.  As  with  the 
previous  test  case,  nonphysical  oscillation  is  a  major  concern.  The  SLl  scheme  generated  much  less  oscilla¬ 
tory  output  than  did  the  NISL  and  Box  schemes.  Diffusion  characteristics  of  Box,  SLl,  and  NISL  are  all 
favorable  compared  to  the  BSBT  scheme.  Figures  16  through  23  show  a  few  comparisons  of  the  BSBT, 
BOX,  NISL,  and  SLl  schemes.  The  SL2  and  SL3  schemes  were  also  tested.  These  schemes  proved  to  be 
less  accurate  than  the  SLl  scheme,  consistent  with  previous  experiments. 

3.3  Computation  Timings 

The  computational  expense  of  the  BSBT,  SLl,  Box,  and  NISL  schemes  were  compared  using  timing 
analysis  for  simple  cases.  The  computational  costs  of  the  four  schemes,  normalized  by  the  BSBT  results, 
are: 


•  BSBT:  1.0 

•  Box:  1.1 

•  SLl:  1.5 

•  NISL:  1.5 

Note  that  for  these  tests,  almost  all  the  computation  time  is  spent  calculating  geographic  propagation. 
By  contrast,  within  SWAN,  most  computation  time  is  spent  solving  the  matrix  required  due  to  finite 
differencing  in  spectral  space.  A  relatively  small  percent  of  time  is  spent  on  the  finite  differencing  of  geo¬ 
graphic  partial  derivatives.  Therefore,  a  50%  increase  in  expense  of  the  geographic  propagation  scheme 
results  in  a  much  smaller  —  approximately  5%  —  increase  in  total  computation  time. 
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Fig.  16  —  Submerged  breakwater  with  gap;  BSBT  seheme.  Wave  aetion  is  speeified  at  the 
left  (x  =  0)  boundary  and  propagated  through  the  domain.  Model  parameters:  CFL  =  0.7; 
direetion  of  transport  theta  =  35  degrees  from  x  axis;  Ax  =  Ay=l  (dimensionless). 
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Fig.  17  —  Submerged  breakwater  with  gap;  Box  seheme.  Model  setup  is  identieal  to 
that  of  the  BSBT  test  ease  deseribed  in  Fig.  16. 
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Fig.  18  —  Submerged  breakwater  with  gap;  NISL  seheme.  Model  setup  is  identieal  to 
that  of  the  BSBT  test  ease  deseribed  in  Fig.  16. 
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Fig.  19 —  Submerged  breakwater  with  gap;  SLl  seheme.  Model  setup  is  identieal 
to  that  of  the  BSBT  test  ease  deseribed  in  Fig.  16. 
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Fig.  20  —  Submerged  breakwater  with  gap;  BSBT  seheme.  Wave  aetion  is  speeified  at 
the  left  (x  =  0)  boundary  and  propagated  through  the  domain.  Model  parameters:  CFL  = 
1.4;  direetion  of  transport  theta  =  35  degrees  from  v  axis;  Ax  =  Ay  =  1  (dimensionless). 
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Fig.  21  —  Test  ease:  submerged  breakwater  with  gap;  Box  seheme.  Model  setup 
is  identieal  to  that  of  the  BSBT  test  ease  deseribed  in  Fig.  20. 
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Fig.  22  —  Test  case:  submerged  breakwater  with  gap;  NISL  scheme.  Model  setup 
is  identical  to  that  of  the  BSBT  test  case  described  in  Fig.  20. 
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Fig.  23  —  Test  case:  submerged  breakwater  with  gap;  SLl  scheme.  Model  setup 
is  identical  to  that  of  the  BSBT  test  case  described  in  Fig.  20. 
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3.4  Conclusions 

The  SLl  scheme  was  chosen  for  implementation  in  SWAN.  The  advantages  of  this  scheme  are: 

•  Numerical  diffusion  is  low  compared  to  the  BSBT  scheme.  SLl  is  also  less  diffusive  than  the  other 
two  Stelling  and  Leendertse  schemes  that  were  tested.  Though  most  of  the  other  numerical  schemes 
tested  exhibited  less  diffusion  than  SLl,  it  is  thought  that  the  diffusion  characteristics  of  SLl  are 
acceptable  for  the  vast  majority  of  applications. 

•  The  scheme  is  unconditionally  stable,  unlike,  e.g.,  the  more  accurate  Backward  Characteristic  Method 
(Lagrangian  Interpolation)  schemes. 

•  The  scheme  has  an  acceptable  computation  cost,  unlike,  e.g.,  the  Borsboom  scheme. 

•  The  scheme  exhibited  relatively  few  problems  with  nonphysical  oscillation,  making  it  favorable 
over  the  Box  scheme  and  several  others. 

4.  IMPLEMENTING  THE  NEW  SCHEME 

The  spatial  propagation  routine  in  SWAN-X33  was  supplemented  with  the  SLl  scheme.  Details  of 
changes  to  the  code  can  be  found  in  Appendix  A.  A  version  of  the  code  that  employs  the  Box  scheme  had 
been  created  previously,  but  development  of  this  code  was  discontinued  when  it  became  clear  that  the  SLl 
scheme  would  be  used  in  SWAN-X33. 

4.1  Scheme  Usage 

Within  three  nodal  points  of  the  geographic  grid  boundaries,  the  BSBT  scheme  is  used.  At  the  interior 
grid  points,  the  user  has  the  option  of  using  either  the  BSBT  scheme  or  the  SLl  scheme. 

The  possibility  of  using  the  Box  scheme  at  the  grid  boundaries  (rather  than  the  BSBT  scheme)  was 
briefly  considered.  However,  since  the  impact  of  this  change  would  be  relatively  minor  (it  would  directly 
influence  only  the  boundary  nodes),  this  alteration  of  the  model  was  not  made. 

4.2  Nonstationary  Mode 

In  nonstationary  mode,  the  sequence  of  operations  in  SWAN-X33  is  identical  to  that  of  SWAN.  The 
BSBT  scheme  is  replaced  by  the  SLl  scheme  at  the  appropriate  grid  locations.  In  one  dimension,  the  simple 
wave  equation  approximated  with  the  SLl  scheme  is 


n-\ 


^•  +  1 


■  -  ^(lOA^;  -  i5A^f_i + 6A^;_  2  -  ni^). 


(45) 


Simple  stepping  in  i  solves  the  matrix.  A  specialized  matrix  solver  is  not  required,  as  all  the  A”  terms 
are  oriented  in  an  upwind  direction,  and  can  therefore  be  treated  independently. 

4.3  Stationary  Mode 

4. 3. 1  Time-Stepping  Method 

Stationary  mode  in  SWAN-X33  cannot  function  as  it  does  in  SWAN  (simple  sweeping),  since  the  SLl 
scheme  is  not  entirely  upwind,  but  is  centered  at  the  explicit  time  level.  Therefore,  the  stationary  model 
time-steps  until  the  simulation  reaches  steady  state  (convergence).  The  user  specifies  the  time  increment,  the 
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maximum  number  of  time  steps,  and  the  eonvergenee  eriteria.  The  user  must  eonsider  the  possibility  of 
nonphysieal  oseillatory  behavior  at  Courant  numbers  signifieantly  greater  than  unity,  and  balanee  this  with 
eonsideration  of  eomputational  eost.  This  is  the  method  eurrently  being  used  by  SWAN-X33.  Alternative 
methods  were  investigated. 

4.3.2  Alternative  Methods 

Two  alternatives  to  the  time-stepping  method  used  in  the  stationary  mode  of  SWAN-X33  were  investi¬ 
gated:  solution  of  the  problem  matrix  by  Gauss-Seidel  method  and  by  relaxation  method.  Setting  the 
terms  in  the  simple  wave  equation  to  zero,  Eq.  (45)  beeomes 

II W,"  =  -  A'"-,')-  ,  +  6NU  -  NU).  (46) 


The  Gauss-Seidel  method  uses  a  slightly  modified  form  of  this  equation — the  time  level  of  the  third  N 
term  is  changed  fi'om  {n-1)  to  («);  since  the  solution  is  a  steady-state  solution,  the  most  current  time  level 
available  is  used  for  all  terms: 


10 
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n: 


-Ci)- 
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(-i5A^;_, + 6aa;_2 


(47) 


Using  the  Gauss-Seidel  method,  the  matrix  solution  did  not  eonverge  for  a  simple,  two-dimensional  test 
ease. 

In  the  relaxation  method,  as  in  the  other  two  methods,  the  matrix  is  solved  by  simple  stepping  in  i.  Using 
the  SLl  seheme  for  the  simple,  one-dimensional  wave  equation,  this  method  is  stated  as 
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(48) 


where  CO  is  a  user-defined  relaxation  parameter.  This  method  was  also  tested  using  a  simple,  two-dimen¬ 
sional  ease.  Using  an  optimal  relaxation  parameter,  eonvergenee  was  aehieved  more  rapidly  than  with  the 
time-stepping  method.  However,  the  eomplexities  of  ehoosing  the  optimal  relaxation  parameter  made  the 
use  of  this  method  in  SWAN-X33  impraetieal. 

5.  VERIFYING  THE  NEW  MODEL 

The  aeeuraey  and  robustness  of  the  model  were  verified  in  a  series  of  tests. 

5.1  Submerged  Breakwater  Test 

As  a  “sanity  eheek”  on  the  new  model,  SWAN-X33  was  eompared  to  the  simple,  two-dimensional  SLl- 
seheme  model  deseribed  in  Seetion  3.2.3.  The  models  were  eompared  in  both  stationary  and  nonstationary 
mode.  The  model  output  was  nearly  identieal  (e.g,.  see  Fig.  24). 
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Fig.  24  —  Test  case:  submerged  breakwater  with  gap;  SLl  scheme;  SWAN-X  is  compared  to  the  simple 
model  used  in  the  scheme  comparisons.  Wave  action  is  specified  at  the  left  (x  =  0)  boundary  and  propagated 
through  the  domain.  Model  parameters:  CFL  =  8e-3;  direction  of  transport  theta  =  30  degrees  from  x  axis; 
Ax  =  Ay  =  1  (dimensionless). 


5.2  Planar  Slope  Test 

Depth-induced  refraction  and  shoaling  in  the  new  model  were  tested  using  a  planar  slope  test  case. 
A  1/125  slope  was  used  with  a  variety  of  wave  conditions.  Results  were  good  (e.g.,  Fig.  25). 

5.3  Tests  with  Currents 

Current-induced  refraction  and  shoaling  were  tested  using  the  same  series  of  experiments  used  by  Ris 
(1997): 

1 .  Following  current:  a  flow  parallel  to  the  wave  direction  increasing  from  0  to  2  m/s. 

2.  Opposing  current:  a  flow  parallel  to  the  wave  direction  increasing  from  0  to  -2  m/s. 

3.  Slanting  current  (0  =  30°):  a  flow  increasing  from  0  to  2  m/s;  wave  angle  of  approach  is  30°,  where 
0  =  -90°  is  a  following  current  (in  the  direction  of  flow)  and  0  =  0  is  cross-flow. 

4.  Slanting  current  (0  =  -30°):  a  flow  increasing  from  0  to  2  m/s;  wave  angle  of  approach  is  -30°. 

Analytical  solutions  for  cases  (1)  and  (2)  are  given  by 


U2  ~  c{c+2U) 


(49) 
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Fig.  25  — A  planar  slope  test  ease.  Model  parameters:  wave  direetion  at  offshore 
boundary  =  30  degrees;  wave  period  =  9.65  s. 


and 


=  i  +  i 
2^2 


(50) 


(e.g.,  Philips  1977).  In  the  slanting  current  case,  if  is  calculated  using  a  conservation  of  energy  relation: 


H^Ho 


sin(20^ ) 
sin(20)  ’ 


(51) 


An  analytical  solution  for  6  was  derived.  Starting  with  the  relations  given  by  Johnson  (1947)  that  as¬ 
sume  that  a  wave  crest  is  continuous  across  a  discontinuity  in  the  current  field,  we  have 


and 


C  ^ 

sin0  sin0 

o 

sin(0  ) 

k  —  k  ^ 

'^~o  sin(0)  ’ 


(52) 


(53) 


and  from  the  dispersion  relation  in  deep  water, 


C  = 


(54) 
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Substituting  Eq.  (53)  into  Eq.  (54)  gives 


C  = 


gsin(0) 
sin(0^ )  • 


(55) 


Substituting  Eq.  (55)  into  Eq.  (52)  gives 


sin  dg 


=  U  + 


sin(e„)sin(e) 


or0  =  arcsin 


1  sin^  d„ 


(56) 


Agreement  with  these  analytieal  solutions  is  good  (see  Figs.  26  through  28). 


Fig.  26  —  Current-induced  refraction  and  shoaling.  Solid  lines  indicate  analytical  solutions; 

“+”  indicates  model  output. 


Results  from  the  tests  of  SWAN  and  SWAN-X  with  eurrent  input  suggest  that  the  models  are  highly 
sensitive  to  frequeney  resolution  (e.g.,  for  the  eases  deseribed  here,  162  frequeney  bins  were  used).  This 
sensitivity  is  thought  to  be  due  to  the  nonzero  Ca  in  the  presenee  of  eurrents.  With  propagation  oeeurring  in 

a-spaee,  the  finite  differeneing  of  the  terms  is  sensitive  to  Ac. 

5.4  Validation  of  Generation  by  Wind 

Wave  generation  by  wind  in  the  modified  model  was  validated  by  eomparison  to  analytieal  solutions 
and  to  the  original  (BSBT-based)  model.  Tests  are  intentionally  similar  to  those  that  were  used  by  Ris  (1997) 
to  validate  the  original  model.  All  simulations  used  stationary,  uniform  wind  fields  of  U^^  =  20  m/s,  identieal 
to  eonditions  used  for  the  SWAMP  (1985)  tests. 
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Fig.  27  —  Current-induced  refraction  and  shoaling;  offshore  wave  angle  =  -30  degrees.  Effect  of 
nonuniform  current  on  wave  direction  is  shown. 


Fig.  28  —  Current-induced  refraction  and  shoaling;  offshore  wave  angle  =  +30  degrees.  Effect  of 
nonuniform  current  on  wave  direction  is  shown. 
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5.4.1  Fetch-Limited  Wave  Growth 

The  model  was  tested  for  eases  of  finite  feteh  lengths,  nominally  infinite  depth,  and  nominally  infinite 
duration.  Nine  simulations,  eaeh  with  a  different  feteh,  were  eondueted:  1124  km,  562  km,  281  km,  141  km, 
70  km,  35  km,  18  km,  8.8  km,  and  4.4  km.  Note  that  the  different  feteh  lengths  eould  have  been  modeled 
using  a  single  simulation,  but  separate  simulations  were  eondueted  to  have  higher  geographie  resolution  for 
the  shorter  fetehes. 

Values  were  nondimensionalized  by  the  frietion  veloeity  as  in  the  SWAMP  tests  of  wind  generation  in 
numerieal  models  (SWAMP  Group  1985): 


Y*  _  8^  -n*  _  8^tot  f*  _  f 

^  ~ur  ~  ut.  g 


(57) 


n2 


where  X*  is  the  dimensionless  feteh;  X  is  the  feteh;  f/*  is  the  frietion  veloeity:  =  -yfc^ 

Cj,  =  (0.8  +  0.065  X  C/jo)  X  0.001  (Wu  1980;  Wu  1982);  is  the  total  energy:  = 

dimensionless  total  energy;  and/^  is  the  dimensionless  peak  frequeney. 


;  is  the 


The  geographie  grid  for  eaeh  of  the  nine  simulations  was  300  x  300  grid  points.  A  relatively  eoarse 
speetral  grid  was  used:  10  frequeney  bins  and  a  direetional  resolution  of  10.3°. 

Results  from  the  two  models  were  eompared  to  empirieal  expressions  formulated  by  Kahma  and  Calkoen 
(1992): 


=  (6.5  X 10“^)  X  (X*)"^  and  4'  =  0.5281  x  (X*)“"'^  (58) 


Agreement  is  fair  (see  Figs.  29  and  30);  it  is  likely  that  agreement  would  have  been  better  if  a  finer 
speetral  resolution  had  been  used  (this  was  not  done  due  to  eomputational  eost). 

Implementation  of  the  new  seheme  seems  to  have  no  deleterious  effeet  on  the  growth  eurve  predietion 
relative  to  that  of  SWAN. 

5.4.2  Depth-Limited  Wave  Growth 

A  seeond  set  of  tests  involved  finite  depth,  nominally  infinite  duration,  and  nominally  infinite  feteh. 
Values  were  nondimensionalized  by  the  veloeity: 


~d-Sd_  e-IE 
ufX-  U!, 


■,andf  = 


4^10 

8 


(59) 


Two  expressions  for  E  and  were  available  for  comparison:  an  analytical  expression  from  Bretschneider 
(1973): 


£  =  5  X 10"'  X 


[tanh(0.53J°’4]",7p  =  0.133  x  [tanh(0.833J°^’4] 


(60) 
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Fig.  29  —  Fetch-limited  wave  generation  by  wind.  Nondimensionalized  energy.  The  BSBT-  and 
SLl -based  models  are  compared  to  an  empirical  growth  curve  from  Kahma  and  Calkoen  (1992). 


Fig.  30  —  Fetch-limited  wave  generation  by  wind.  Nondimensionalized  frequency.  The  BSBT-  and 
SLl -based  models  are  compared  to  an  empirical  growth  curve  from  Kahma  and  Calkoen  (1992). 
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and  an  empirical  expression  by  Young  and  Verhagen  (1996),  which,  given  infinite  fetch,  is  written  as: 

E  =  3.64  X  lO"'  X  [tanh(0.493J°’')]‘  '",4  =  0-133  x  [tanh(0.331J“’‘)]"°'\  (61) 

The  latter  set  of  equations  was  derived  for  (6  X 10”^)  <d  <2. 

Seven  simulations  were  conducted  (see  Table  1).  The  time  step  used  in  SLl  computations  is  also  shown. 
Note  that  the  Courant  numbers  for  these  simulations  were  high  [0(100)].  Large  time  steps  were  required  to 
approximate  infinite  duration  while  keeping  computational  time  reasonable.  Due  to  the  uniform  and  steady 
model  input  (both  wind  field  and  bathymetry),  the  oscillations  in  the  wave  field  (resulting  from  high  Cou¬ 
rant  number)  were  not  severe. 


Table  1  —  Simulation  Parameters  for  Depth-Limited  Wave  Growth 


Depth  (m) 

Ax  (m) 

At  (s)  (SLl  only) 

Approximate  q  (m/s) 

0.45 

10 

2000 

1 

1.33 

30 

2000 

2 

4 

60 

2000 

4 

12 

375 

3000 

7 

36 

1250 

15000 

10 

108 

2500 

25000 

11 

324 

3750 

40000 

11.5 

The  Hasselmann  (1973)  formulation  for  bottom  friction,  with  a  JONSWAP  friction  coefficient  of 
C^^^=0.038m^s'^  was  used  for  all  simulations.  Hasselmann  (1973)  suggests  a  of  0.038m^s'^  for  swell 
conditions  and  Bouws  and  Komen  (1983)  propose  a  of  0.067m^s'^  for  wind  sea  conditions.  For  these 
simulations,  was  not  varied,  as  Ris  (1997)  observed  that  these  tests  are  not  sensitive  to  friction  coeffi¬ 
cient.  Comparisons  of  both  models  to  the  analytical  and  empirical  solutions  are  good  (see  Figs.  31  and  32). 

5.4.3  Duration-Limited  Wave  Growth 

For  duration-limited  wave  growth,  both  infinite  depth  and  infinite  fetch  are  approximated  for  several 
simulation  durations:  83  min,  3  hr,  14  hr,  28  hr,  56  hr,  and  139  hr.  The  two  models  are  compared  to  verify 
that  the  higher-order  scheme  does  not  have  any  impact  on  growth  rate.  Indeed,  this  is  the  case  (see  Fig.  33). 

5.5  Berkhoff-Booij-Radder  (BBR)  Shoal 

Berkhoff  et  al.  (1982)  conducted  a  series  of  experiments  of  wave  propagation  over  a  submerged  ellipti¬ 
cal  shoal  superimposed  on  a  plane  beach.  This  experimental  setup  provides  a  scenario  where  both  refraction 
and  diffraction  are  significant. 
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Fig.  3 1  — Fetch-limited  wave  generation  by  wind.  Nondimensionalized  energy.  The  BSBT- 
and  SLl -based  models  are  compared  to  growth  curves  taken  from  Bretschneider  (1973)  and 
Young  and  Verhagen  (1996). 
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Fig.  32  —  Fetch-limited  wave  generation  by  wind.  Nondimensionalized  frequency.  The 
BSBT-  and  SLl -based  models  are  compared  to  growth  curves  taken  from  Bretschneider 
(1973)  and  Young  and  Verhagen  (1996). 
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Fig.  33  —  Duration-limited  wave  generation  by  wind.  Wave  heights  are  plotted. 
The  BSBT-  and  SLl -based  models  are  eompared. 


Domain  size:  25  x  25  m 
Grid  resolution:  Nx  =  Ay  =  0.25  m 
Wave  height  at  boundary:  =  0.0464  m 

Wave  period:  7=  1  s  (monoehromatie) 

Wave  angle  at  input  boundary:  0  =  0°  (unidireetional) 

Figure  34  shows  the  experimental  layout.  Data  from  the  experiment  were  eompared  to  results  from 
SWAN,  SWAN-X33,  and  REF/DIFl  (Kirby  and  Dalrymple  1993)  —  a  phase-resolving,  monoehromatie 
wave  model  based  on  the  parabolie  approximation  of  the  mild-slope  equation.  Unlike  the  two  SWAN  mod¬ 
els,  REF/DIFl  has  no  diffusion  and  ineludes  the  effeets  of  both  refraetion  and  diffraetion  (implieitly).  In 
these  simulations,  refraetion  ereates  a  feature  in  the  wave  height  field  that  is  prone  to  diffusion  in  the  hyper- 
bolie  models.  SWAN-X33  results  show  signifieantly  less  diffusion  than  SWAN.  Figures  35  and  36  show 
results  along  two  transeets  where  the  differenee  between  the  two  models  is  greatest  (transeet  nos.  3  and  7). 
Note  the  absenee  of  diffraetion-indueed  oseillations  in  H(y)  output  from  the  SWAN  models. 

5.6  DELILAH 

5. 6. 1  Introduction 

The  DELILAH  (Duek  Experiment  on  Low-frequeney  and  Ineident-band  Longshore  and  Aeross-shore 
Hydrodynamies)  (Birkemeier  1991)  was  eondueted  at  the  U.S.  Army  Corps  of  Engineers  Field  Researeh 
Faeility  at  Duek,  North  Carolina  during  Oetober,  1990  in  order  to  investigate  many  aspeets  of  surf  zone 
physies.  Numerous  measurement  instruments  were  employed  and  several  researeh  entities  partieipated. 
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Fig.  34  —  Layout  of  the  Berkhoff-Booij-Radder  lab  experiment 
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Fig.  35  —  The  BBR  test  ease:  eomparison  of  data  to  the  phase-resolving  model 
REF/DIFl  and  the  two  SWAN  models  along  transeet  3  (a  lateral  transeet) 


Improving  the  Numerics  of  a  Third-Generation  Wave  Action  Model 


41 


2.5 


2- 


1.5- 


o 

X 


1 


0.5- 


p  o  o 


o  measured 

.  REF/DIF1  linear 

-  -  SWAN 
■  -  ■  -  SWAN-X  w/SLI 


10 


12 


14 


16 

x(m) 


18 


20 


22 


Fig.  36  —  The  BBR  test  ease:  eomparison  of  data  to  the  phase-resolving  model 
REF/DIFl  and  the  two  SWAN  models  along  transeet  7  (a  longitudinal  transeet) 


5. 6.2  Available  Data 

Daily  surveys  of  a  550  m  by  375  m  minigrid  area  were  eombined  with  a  single,  larger  (1590  m  x  1640 
m)  regional  grid,  providing  model  bathymetries  for  eaeh  day  of  the  experiment.  Bathymetrie  eontours  are 
approximately  straight  and  parallel,  exeept  near  a  pier  (whieh  is  outside  the  minigrid  area).  At  Duek,  sand 
bar  features  are  often  formed  near  the  shoreline  during  storm  events. 

Nearshore  wave  gauges  and  eurrent  meters  eolleeted  data  inside  the  minigrid  area.  The  nearshore  wave 
gauges  provide  direetional  speetra.  However,  beeause  of  less  than  optimal  plaeement  of  the  wave  gauges, 
the  nearshore  speetra  are,  for  praetieal  purposes,  nondireetional.  The  nearshore  gauges  also  provide  waveheight 
data  along  one  eross-shore  transeet  and  two  longshore  transeets. 

High-quality  direetional  speetra  data  were  available  for  the  permanent  8-m  depth  “Linear  Array”  and 
for  the  13-m  depth  SAMSON  array.  The  models  were  initialized  using  two-dimensional  speetra  that  had 
been  ealeulated  from  the  SAMSON  array  measurements.  Data  were  taken  at  2  Hz  and  samples  were  taken 
for  2  hr,  16  min.  periods.  These  samples  represented  3  hr  time  spans,  with  the  remaining  44  min.  taken  for 
data  analysis  and  proeessing. 

Figure  37  shows  the  wave  gauge  arrangements  and  a  general  map  of  the  experiment  area. 

5.6.3  Model  Results 

SWAN  and  SWAN-X33  were  eaeh  run  for  15  of  the  3-hr  periods  between  Oetober  6  and  27,  1990,  to 
test  a  full  range  of  wave  eonditions.  The  bathymetry  itself  also  ehanged  eonsiderably  during  this  time.  In 
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Fig.  37  —  Layout  of  the  DELILAH  simulations.  Locations  of  spectral  data  are  shown  as  x’s. 
Locations  of  nearshore  gauges  are  denoted  with  o’s.  Depth  contours  are  in  meters.  NNW  is  up. 


both  models,  nonlinear  wave-wave  interaetions  were  turned  off  and  wave  generation  by  wind  was  not  in- 
eluded.  Dissipation  by  bottom  frietion  was  also  omitted  (this  is  the  default  setting  of  the  model). 

These  15  eases  were  also  modeled  using  the  speetral  parabolie  wave  model  REF/DIF-S  (Kirby  and 
Ozkan  1994)  and  a  simple  model  that  ealeulates  wave  transformation  using  the  approximation  of  straight 
and  parallel  eontours.  In  general,  all  four  wave  models  provided  very  similar  results:  model-to-model  differ- 
enees  were  typieally  mueh  smaller  than  model-to-data  differenees. 

5.6.3.1  Global  Waveheights 

The  differenee  between  SWAN  and  SWAN-X  eaused  by  diffusion  is  most  easily  visualized  by  eompar- 
ing  global  waveheights.  Refraetion-indueed  features  are  evident,  espeeially  those  ereated  by  the  seour  hole 
in  the  vieinity  of  the  Field  Researeh  Faeility  pier  (see  Fig.  38).  These  features  are  visibly  diffused  by  the 
BSBT  seheme.  Note  that  even  without  diffusion,  these  refraetion  features  are  fairly  smooth  due  to  the  aver¬ 
aging  effeet  of  a  full  wave  speetrum.  Spatial  gradients  and  the  resulting  diffusion  are  both  redueed  by  the 
broadness  of  the  input  speetra.  For  more  unidireetional  and  monoehromatie  wave  eonditions,  the  impaet  of 
the  numerieal  diffusion  would  be  mueh  greater. 

5.6.3.2  Waveheights  at  Cross-Shore  Transects 

Beeause  diffusion  tends  to  aet  perpendieularly  to  the  axis  of  wave  when  eomputations  are  stationary, 
eomparisons  of  the  two  SWAN  models  along  a  eross-shore  transeet  do  not  effeetively  illustrate  diffusion. 
However,  sinee  high-quality  data  were  available  along  sueh  a  transeet,  eomparisons  were  made  as  a  valida¬ 
tion  of  the  two  models.  Agreement  with  data  is  generally  good  for  the  models  (e.g.,  see  Fig.  39),  though 
SWAN  overprediets  waveheights  in  the  nearshore  for  most  of  the  eases.  This  ean  be  partially  attributed  to 
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Fig.  38  —  Global  waveheight  comparison  for  DELILAH  simulation  of  10-07-90@1600  (EST). 
Zero  moment  waveheights  are  shown  in  meters. 


Fig.  39  —  Zero  moment  waveheights  are  compared  along  a  cross-shore  transect: 
DELILAH  simulation  of  10-14-90@1900  (EST) 
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the  omission  of  bottom  friction,  although  other  factors  may  have  been  equally  important  (e.g.,  possible 
deficiencies  of  SWAN’s  breaking  mechanism). 

5.6.3.3  Directional  Spectra  at  8-m  Contour 

Model  output  was  compared  to  data  from  the  8-m  gauge  location.  As  one  might  expect,  the  degree  of 
wave  transformation  between  the  13-m  depth  contour  and  8-m  depth  contour  is  generally  quite  small. 
Figures  40  through  43  show  the  spectral  comparisons  for  two  simulations  with  a  typical  error  level  (error  in 
terms  of  measured  vs  computed).  Both  models  compare  well  with  the  data.  SWAN  and  SWAN-X33 
show  very  similar  results,  except  at  one  end  of  the  directional  spectra,  where  the  SWAN-X  spectra  are 
somewhat  irregular.  This  irregularity  is  typical  of  the  15  SWAN-X  simulations  conducted.  The  irregularity 
is  caused  by  interaction  of  the  SLl  scheme  in  geographic  space  with  the  centered  difference  scheme  in 
directional  space.  It  is  thought  that  this  irregularity  has  minimal  deleterious  impact  on  model  results,  par¬ 
ticularly  in  regards  to  quantities  like  waveheights. 


Fig.  40  —  Comparison  of  frequency  spectra  for  the  two  models  and  data  at  the  8  m 
array  for  DELILAH  simulation  of  10-10-90@0100  (EST) 


5.6.3.4  Nondirectional  Spectra  at  Nearshore  Gauges 

Measured  frequency  spectra  from  the  two  nearshore  arrays  were  compared  to  model  results.  Figures  44 
and  45  show  typical  results  for  the  (seaward)  bar  crest  and  the  (landward)  bar  trough  array,  respectively. 
Comparisons  of  the  SWAN  models — BSBT  mode  vs  SLl  mode — show  an  expected  outcome:  similar  S(f) 
distributions  with  small  differences  in  waveheight  due  to  decreased  diffusion  of  refraction-focusing  in  the 
SLl -based  model. 
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Fig.  41  —  Comparison  of  directional  spectra  for  the  two  models  and  data  at  the  8-m 
array  for  DELILAH  simulation  of  10-10-90@0100  (EST) 


Fig.  42  —  Comparison  of  frequency  spectra  at  the  8-m  array  for  DELILAH  simulation 
of  10-11-90@0100  (EST):  data  vs  SWAN,  REF/DIF-S,  and  simplified  model 
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Fig.  43  —  Comparison  of  directional  spectra  at  the  8-m  array  for  DELILAH  simulation 
of  10-11-90@0100  (EST):  data  vs  SWAN,  REF/DIF-S,  and  simplified  model 


Fig.  44  —  Comparison  of  frequency  spectra  at  the  bar  crest  array  for  DELILAH  simulation 
of  10-14-90@1900  (EST):  data  vs  the  two  SWAN  models 
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Fig.  45  —  Comparison  of  frequency  spectra  at  the  bar  trough  array  for  DELILAH 
simulation  of  10-14-90(gl900  (EST):  data  vs  the  two  SWAN  models 


Although  not  particularly  relevant  to  a  discussion  of  numerics  and  diffusion,  the  effect  of  triad  wave- 
wave  interaetion  on  model  results  is  worth  mentioning.  Beeause  amplifieations  of  the  seeond  harmonie  of 
the  speetral  peak  were  observed  in  the  nearshore  data,  several  DELILAH  eases  were  simulated  with  triads 
enabled  in  SWAN.  Generally,  this  resulted  in  a  large  improvement  in  model-data  eomparisons  (e.g.,  see  Lig. 
46).  However,  the  model  still  overprediets  wave  energy  in  the  nearshore,  partieularly  at  the  higher  frequen- 
eies. 

5.7  DUCK94 

5. 7. 1  Introduction 

The  DUCK  94  experiment  (Birkemeier  and  Thornton  1994)  took  plaee  in  August  and  Oetober  of  1994. 
The  primary  purpose  was  to  study  sediment  transport,  morphology,  wave  transformation,  nearshore  eireula- 
tion,  and  swash  proeesses,  as  well  as  to  provide  a  meehanism  for  eollaboration  between  researehers.  Nine¬ 
teen  organizations  and  over  100  researehers  partieipated.  Wave  data  eolleeted  as  part  of  the  DUCK94  ex¬ 
periment  were  used  to  validate  and  eompare  the  two  SWAN  models. 

5. 7.2  Model  Setup  and  Input 

The  models  were  initialized  with  data  eolleeted  by  the  “inner  shelf  buoy”  that  had  been  deployed  for  the 
experiment  (Jensen  1995).  Proeessed  wave  speetra  at  this  buoy  were  available  hourly  for  the  time  period  of 
0000  GMT  Get.  10  to  0000  GMT  Get.  22.  Wind  data,  also  eolleeted  at  this  buoy,  were  used  for  some  of  the 
simulations.  The  buoy  was  loeated  in  approximately  27-m  water  depth. 
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Fig.  46  —  Comparison  of  frequency  spectra  at  the  bar  crest  array  for  DELILAH  simulation 
of  10-14-90@1900  (EST):  data  vs  SWAN  (with  and  without  triads) 


Spectral  output  from  the  models  was  compared  to  data  collected  at  the  U.S.  Army  Corps  of  Engineers 
Field  Research  Facility  Linear  Array,  located  in  approximately  8-m  water  depth.  Directional  wave  spectra 
are  calculated  from  data  collected  at  this  permanent  array  every  three  hours.  During  the  12  days  that  were 
modeled,  there  are  89  time  periods  for  which  both  inner  shelf  data  and  Linear  Array  data  exist.  Seven 
simulation  sets  (each  composed  of  89  simulations)  were  conducted.  Four  are  relevant  to  this  study: 

•  SWAN-X  with  the  BSBT  scheme,  without  wind 

•  SWAN-X  with  the  SLl  scheme,  without  wind 

•  SWAN-X  with  the  BSBT  scheme,  with  wind 

•  SWAN-X  with  the  SLl  scheme,  with  wind 

These  simulations  did  not  include  bottom  friction. 

The  Coastal  Hydraulics  Laboratory,  U.S.  Army  Corps  of  Engineers,  provided  a  bathymetric  grid  for  the 
region.  This  grid  was  created  by  Dr.  John  Dugan  (of  Arete  Associates)  using  NOS  digital  bathymetry,  DBDB5 
bathymetry,  and  hydrographic  surveys  conducted  during  DUCK94.  The  grid  has  the  following  characteris¬ 
tics: 

•  40.0  km  across  in  x  (longitude). 

•  107  km  wide  in  y  (latitude). 

•  404  meter  resolution  in  x. 

•  500  meter  resolution  in  y. 

•  Low-x  side  of  the  grid  corresponds  to  the  longitude  of  the  Inner  Shelf  buoy;  wave  height  is  specified 
along  this  boundary. 

•  High-x  side  of  the  grid  is  land. 

Figure  47  shows  the  bathymetry  and  the  locations  of  the  instruments. 
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Fig.  47  —  Bathymetry  used  for  DUCK94  simulations.  Depths  are  shown  in  meters. 


For  the  simulations  with  the  higher  order  seheme,  speeifieation  of  a  time  step  size  was  neeessary.  This 
value  was  ehosen  sueh  that  the  peak  wave  of  a  partieular  speetrum  would  propagate  aeross  the  domain  four 
times  during  the  simulation.  The  approximate  travel  time  of  the  peak  wave  was  approximated  using  the 
direetion  and  group  veloeity  of  the  wave  over  an  approximated  flat  bathymetry.  The  ealeulated  time  step 
sizes  resulted  in  CFLs  ranging  from  8  to  15. 

5.7.3  Comparison:  Simulations  Without  Wind  Input 

Differenees  in  global  waveheights  produeed  by  the  two  models  are,  for  the  most  part,  unremarkable. 
This  is  not  surprising,  as  model  input  is  stationary  and  spatial  variability  of  wave  eonditions  is  minimal. 
Refraetion-indueed  features  do  exist,  however,  and  diffusion  of  these  features  is  apparent.  Figure  48  shows 
waveheights  in  the  region  of  the  Linear  Array.  The  more  distinetive  refraetion  features  in  the  SLl -based 
model  suggest  that  the  BSBT  seheme  is  artifieially  diffusing  these  features. 

The  impaet  of  the  numerieal  seheme  on  wave  eonditions  at  the  loeation  of  the  Linear  Array  is  minimal. 
Figures  49  through  52  eompare  gross  wave  parameters  for  both  models  and  data.  The  two  models  show 
similar  results;  both  eompare  fairly  well  to  the  data,  though  admittedly,  these  simulations  of  wave  propaga¬ 
tion  are  not  diffieult  tests  of  model  skill. 

Direetional  speetra  at  the  Linear  Array  (e.g.,  as  in  Figs.  53  and  54)  were  eompared,  but  no  elear  trends 
were  observed.  Error  is,  in  general,  evenly  distributed  aeross  the  speetra. 
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Fig.  48  —  Comparison  of  global  waveheights  from  the  two  SWAN  models  on  10-19-1994@0600  EST. 
Waveheights  have  been  normalized  by  the  waveheights  at  the  Inner  Shelf  buoy. 
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Fig.  49  —  Comparison  of  waveheights  from  the  two  SWAN  models  at  the  linear 
array  for  the  12  days  that  were  modeled  (zero  moment  waveheights  in  meters) 
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Fig.  50  —  Comparison  of  normalized  waveheights  from  the  two  SWAN  models  at  the  linear 
array  for  the  12  days  that  were  modeled.  Waveheights  have  been  normalized  by  the  waveheights 
at  the  Inner  Shelf  buoy. 


Fig.  5 1 - Comparison  of  peak  wave  angles  from  the  two  SWAN  models  at  the  linear  array 

for  the  12  days  that  were  modeled.  Wave  angles  are  in  degrees  (traveling  toward  north  =  -90 
degrees;  traveling  toward  south  =  +90  degrees) 
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Fig.  52  —  Comparison  of  peak  wave  periods  (in  seeonds)  from  the  two  SWAN 
models  at  the  linear  array  for  the  12  days  that  were  modeled 


Fig.  53  —  Comparison  of  frequeney  speetra  at  the  8-m  array:  data  vs  the 
two  SWAN  models.  Simulation:  10-14-1994@1800  EST. 
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direction  (deg) 


Fig.  54  —  Comparison  of  directional  spectra  at  the  8-m  array:  data  vs  the 
two  SWAN  models.  Simulation:  10-14-1994@1800  EST. 


5.7.4  Comparison:  Simulations  With  Wind  Input 

Comparisons  of  gross  wave  parameters  to  the  data  are  visibly  improved  with  the  inelusion  of  wind  for 
both  the  SWAN  simulations  (Fig.  55)  and  the  SWAN-X  with  SLl  simulations  (Fig.  56).  Unsurprisingly,  the 
improvement  is  most  notieeable  during  the  first  two  days  of  the  data  set,  when  loeal  wind  sea  dominates.  As 
in  the  simulations  without  wind,  minor  differenees  between  the  two  sehemes  exist  in  the  global  waveheights, 
but  at  the  linear  array,  the  impaet  of  the  redueed  diffusion  is  insignifieant. 

5.8  Southern  California  Bight 

5. 8. 1  Introduction 

San  Miguel  and  Santa  Rosa  Islands  are  loeated  offshore  of  Ventura,  California,  in  the  Southern  Califor¬ 
nia  Bight.  Figure  57  shows  their  loeation  with  respeet  to  the  eoast.  The  wave  elimate  in  the  Southern  Califor¬ 
nia  Bight  (stretehing  from  Point  Coneeption  south  to  the  Mexiean  border)  has  been  aetively  monitored  by 
the  Coastal  Data  Information  Program  (CDIP),  whieh  is  based  at  Seripps  Institution  of  Oeeanography  in  La 
Jolla.  CDIP  is  jointly  funded  by  the  California  Department  of  Boating  and  Waterways  and  the  U.S.  Army 
Corps  of  Engineers. 

The  California  shelf  is  extremely  narrow;  it  is  generally  no  more  than  1 1  km  from  shoreline  to  shelfbreak. 
The  bathymetry  in  the  Bight  is  mostly  planar,  with  submarine  eanyons  and  other  features  in  some  areas.  The 
California  eoastline  south  of  Point  Coneeption  generally  faees  southwest.  The  presenee  of  Point  Coneeption 
shelters  mueh  of  the  Bight  from  waves  eoming  from  the  northwest,  as  oeeurs  often  during  winter  storms. 
During  the  summer,  the  wave  elimate  in  the  Bight  is  dominated  by  swell  generated  by  Southern  Hemisphere 
storms.  The  islands  in  the  Bight  potentially  shelter  mueh  of  the  eoastline  from  offshore  waves. 
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Fig.  55  —  Comparison  of  normalized  waveheights  from  SWAN-X  with  BSBT  scheme  (with 
and  without  wind)  at  the  linear  array  for  the  12  days  that  were  modeled.  Waveheights  have 
been  normalized  by  the  waveheights  at  the  Inner  Shelf  buoy. 
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Fig.  56  —  Comparison  of  normalized  waveheights  from  SWAN-X  with  SLl  seheme  (with 
and  without  wind)  at  the  linear  array  for  the  12  days  that  were  modeled.  Waveheights  have 
been  normalized  by  the  waveheights  at  the  Inner  Shelf  buoy. 
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West  Longitude 


Fig.  57  —  Map  of  San  Miguel  (left)  and  Santa  Rosa  (right)  Islands  and  vieinity 


In  early  1992,  pressure  gauges  were  installed  around  Santa  Rosa  Island  to  measure  the  wave  elimate  in 
the  area.  Figure  58  shows  the  loeations  of  the  gauges  and  the  bathymetry  around  the  islands.  We  will  empha¬ 
size  eomparisons  of  the  data  with  the  model  at  gauge  #10  (loeated  between  the  islands)  and  gauge  #11 
(loeated  on  the  southwest  side  of  Santa  Rosa  Island).  A  Datawell  buoy  loeated  offshore  of  Point  Coneeption 
measured  waves  as  they  entered  the  Bight. 

5.8.2  Model  Setup 

We  used  the  bathymetry  depieted  in  Fig.  58  to  simulate  swell  propagation  around  the  islands.  The 
presenee  of  the  islands  would  tend  to  exaeerbate  the  diffusion  in  the  geographie  propagation  terms  in  SWAN. 
The  bathymetrie  grid  extended  66.7  km  by  53.4  km,  with  a  resolution  of  100  meters  in  both  direetions.  This 
very  fine  resolution  was  used  to  allow  direet  eomparisons  to  a  phase-resolving  model;  it  ensured  that  the 
bathymetrie  variations  in  the  vieinity  of  the  islands  would  be  well  resolved  in  the  model  (numerieal  experi¬ 
ments  later  eonfirmed  a  signifieant  dependenee  on  resolution  at  some  loeations).  Direetional  speetra  mea¬ 
sured  at  the  Point  Coneeption  buoy  during  January  1992  were  used  as  model  input,  with  the  frequeneies 
limited  to  the  swell  range  (0.05  to  0.09  Hz).  The  direetional  range  used  eomprised  waves  arriving  from  due 
north  to  due  south  (37  direetions  at  5  degree  intervals). 

5.8.3  Model  Results 

The  stationary  SWAN  model  was  run  with  all  January  1992  speetra.  Figure  59  eompares  the  energy  (in 
em^)  from  SWAN  output  with  that  of  the  gauge  data  at  gauge  #10.  Also  shown  are  results  from  a  speetral 
refraetion  model  (LeMehaute  and  Wang  1982)  and  a  parabolie  refraetion/diffraetion  model  (Kirby  1986). 
The  implementation  of  the  latter  two  models  for  the  Southern  California  Bight  was  detailed  in  O’Reilly  and 
Guza  (1993).  Finally,  the  energy  level  at  the  Point  Coneeption  buoy  is  also  shown.  The  inelusion  of  the 
offshore  data  in  the  plot  emphasizes  the  sheltering  experieneed  by  gauge  #10.  In  general,  SWAN  overprediets 
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Fig.  58  —  Bathymetry  and  gauge  loeations  around  San  Miguel  and  Santa  Rosa  Islands. 
Depths  are  in  meters.  Gauges  shown  are:  #8  (o),  #9  (x),  #10  (*),  and  #11  (+). 


Day  in  January  1992 


Fig.  59  —  Comparisons  of  SWAN  energy  predietions  (*)  at  gauge  #10  to  measured  data  (o), 
speetral  refraetion  model  (x),  and  refraetion/diffraetion  model  (+).  Dashed  line  is  energy 
level  at  Point  Coneeption  buoy. 
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the  energy  at  this  gauge.  This  is  a  eonsequenee  of  the  diffusive  numerieal  seheme  leaking  energy  into  heavily 
sheltered  areas.  Figure  60  shows  a  eomparison  of  the  models  to  the  measurements  at  gauge  #11.  This  site  is 
mueh  more  exposed  to  the  offshore  wave  elimate.  For  the  most  part,  the  SWAN  model  does  as  well  as  the 
other  models  in  estimating  the  wave  energy. 


Fig.  60  —  Comparison  of  SWAN  energy  predietions  (*)  at  gauge  #1 1  to  measured  data  (o), 
speetral  refraetion  model  (x),  and  refraetion/diffraetion  model  (+) 


Due  to  the  fine  geographie  resolution  used  (with  668  x  535  grid  points),  the  simulations  were  expensive 
eomputationally.  The  SWAN  runs  for  the  entire  January  1992  data  set  required  a  total  of  approximately 
2  days  of  CPU  time  on  a  Sun  Ultra2300  workstation  (300  MHz  proeessor).  Eaeh  run  required  approximately 
300  MB  of  RAM.  Beeause  of  the  fine  spatial  resolution,  an  equally  fine  temporal  resolution  was  needed  due 
to  CFL  eonsiderations:  the  SWAN-X  simulations  with  the  SLl  seheme  used  a  At  of  10  seeonds.  With  this 
small  time  step  inerement,  1000  time  steps  were  neeessary  in  order  to  propagate  energy  to  a  steady-state 
eondition.  Thus,  each  SWAN-X  run  required  almost  2  weeks  of  CPU  time  and  approximately  600  MB 
RAM  on  the  same  workstation.  Clearly,  use  of  the  pseudo-stationary  SLl  model  with  a  large  (100,000+  grid 
points),  high-resolution  (less  than  200  m  spaeing)  grid  is  impraetieal  for  any  present  workstation  or  PC. 
Fortunately,  unlike  phase-resolving  models,  SWAN  does  not  require  a  minimum  spatial  resolution,  so  when 
the  bathymetry  does  not  require  high  resolution,  the  size  of  the  SWAN-X  simulation  ean  be  tailored  to 
eomputational  resourees. 

The  results  do  look  promising  in  the  vieinity  of  gauge  #10.  Figure  61  shows  frequeney  and  direetion 
speetra  at  gauge  #10  from  a  SWAN  run  for  January  14,  1992  at  2300  PST.  The  direetion  speetra  S(Q)  shows 
two  peaks,  eorresponding  to  energy  propagating  around  both  sides  of  San  Miguel  Island.  (Zero  degrees 
indieates  due  east  propagation.)  The  zero  moment  waveheight  =  1.03  m.  Figure  62  shows  the  same 
result  from  SWAN-X.  The  zero  moment  waveheight  is  notieeably  lower:  0.92  m.  Table  2  summarizes  the 
results  for  this  ease  at  gauge  #10. 
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SWAN  Hmo=1.03m  Tp=1 7.241 4s  thetap=35deg 


Fig.  61  — Frequency  (top)  and  direction  (bottom)  spectra  from  SWAN  at  gauge  #10  for 
January  14,  1992  at  2300  PST 


SWAN-X  Hmo=0.92m  Tp=1 7.241 4s  thetap=40deg 


Fig.  62  —  Frequency  (top)  and  direction  (bottom)  spectra  from  SWAN-X  at  gauge 
#10  for  January  14,  1992  at  2300  PST 
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Table  2  —  Comparisons  of  at  Gauge  #10  for  January  14,  1992 

at  2300  PST 


at  Gauge  #10  (m) 

Measured 

0.89  m 

Refraction/Diffraction 

0.94  m 

Refraction 

0.92  m 

SWAN 

1.03  m 

SWAN-X 

0.92  m 

The  reduction  in  the  SWAN-X  H  prediction  relative  to  that  of  SWAN  is  due  to  the  reduced  numerical 
diffusion  of  the  higher  order  seheme.  Globally,  the  effeet  of  this  seheme  ean  be  seen  from  the  waveheight 
predietions  throughout  the  domain  (Fig.  63).  The  SWAN  predietion  shows  a  fairly  smooth  waveheight  field 
on  the  northern  side  of  Santa  Rosa  Island  and  between  it  and  San  Miguel  Island.  This  smoothness  is  an 
artifaet  of  numerieal  diffusion.  In  eontrast,  the  SWAN-X  predietion  shows  mueh  more  variability  in  the 
waveheight  patterns  around  the  island,  due  to  the  higher-order  seheme  for  geographieal  propagation. 

6.  THE  GARDEN  SPRINKLER  EFFECT 

The  “garden  sprinkler  effeet”  (GSE)  is  a  well-known  problem  whereby  an  initial  spatial  distribution  of 
wave  energy  disintegrates  into  individual  eomponents  as  the  energy  propagates,  eaeh  eomponent  eorre- 
sponding  to  a  partieular  speetral  bin  used  in  the  model.  The  GSE  is  perhaps  most  easily  deseribed  as  a 
problem  of  speetral  resolution:  with  an  extremely  fine  speetral  resolution,  energy  would  disperse  evenly  and 
the  resulting  wave  field  would  be  —  in  most  instanees  —  eontinuous.  With  a  eoarse  speetral  resolution,  the 
resulting  energy  field  is  diseontinuous  due  to  geographie  separation  of  energy  bins  (e.g.,  see  Figs.  64  and 
65).  Typieally,  eomputational  resourees  require  the  use  of  a  relatively  eoarse  speetral  resolution  (A6  >  5°). 
Beeause  the  problem  is  greater  for  eases  of  severe  spatial  and  temporal  variability  of  input,  it  is  generally 
assoeiated  with  propagation  at  large  seales. 

Numerieal  diffusion  tends  to  eounter  the  GSE,  as  it  smooths  the  spatial  irregularities  in  the  energy  field. 
Thus,  the  problem  is  most  notieeable  in  models  that  have  little  numerieal  diffusion.  With  the  introduetion  of 
the  new  propagation  seheme  in  SWAN-X,  there  is  a  greater  ineentive  for  ineorporating  measures  to  eounter 
the  GSE. 

Booij  and  Holthuijsen  (1987)  ereated  one  sueh  tool.  It  operates  by  diffusing  wave  energy  a)  in  the 
direetion  of  wave  propagation  (to  eounter  effeets  of  frequeney  diseretization)  and  b)  normal  to  the  direetion 
of  wave  propagation  (to  eounter  effeets  of  direetional  diseretization).  The  magnitude  of  the  diffusion  for 
eaeh  is  eontrolled  by  a  diffusion  eoeffieient — D  and  D  ,  respeetively.  The  eoeffieients  are  a  fiinetion  of  the 
level  of  diseretization  and  the  propagation  eharaeteristies  of  the  wave  energy: 

A,=Acfr/12 
D„„=cfA0fr/12 


(62) 
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Hmo  (m):  SWAN-X  w/BSBT  scheme 


Hmo  (m):  SWAN-X  w/SL1  scheme 
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Fig.  63  —  Global  waveheight  results  from  SWAN  (top)  and  SWAN-X  (bottom)  for 
January  14,  1992  at  2300  PST 


X  (dimensionless) 


Fig.  64  —  An  energy  field.  All  energy  starts  at  the  same  loeation,  shown  in 
lower  left  eomer  (x=y  =  3).  Energy  disperses  as  it  propagates  due  to  direetional 
and  frequeney  (i.e.,  speed)  distribution  of  the  initial  speetra.  In  this  ease,  the 
initial  speetrum  is  very  finely  resolved,  so  the  Garden  Sprinkler  Effeet  is  insig- 
nifieant.  The  resulting  energy  field  is  eontinuous. 
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Fig.  65  — An  energy  field.  The  diseontinuity  is  a  result  of  the  Garden  Sprinkler  Effeet. 

All  energy  starts  at  the  same  loeation,  shown  in  lower  left  eomer  (x  =  y  =  3).  Energy 
disperses  as  it  propagates  due  to  direetional  and  frequeney  (i.e.,  speed)  distribution  of 
the  initial  speetrum.  In  this  ease,  the  initial  speetrum  eonsists  of  8  frequeney  bins  and  8 
direetional  bins. 


In  polar  coordinates,  the  anti-GSE  equation  is  written  as: 

dt  as  ds  ds  dndn 

Converted  to  Cartesian  coordinates  (in  which  SWAN  is  written), 


3  3  3  3  3 

—  N^■  {x,  +  ^  {x,  y,  0  -  ^  ^  0]  +  ^  [c^yNj  y’ 

-2Z)^|-y,(x,y,0  =  0, 

where 

A.  =  A,  cos' 0  + A.  sin- e 
=  D„  sin' e  +  D„.  cos' 0 
A,  =  (A,-A.)sin0cos0 


(65) 
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Here,  “s'"  denotes  frequeney  diseretization,  “n”  denotes  direetional  diseretization,  N  is  the  average 
aetion  density  in  a  partieular  aetion  bin,  and  T  is  “wave  age,”  or  the  time  elapsed  sinee  generation/initializa¬ 
tion  of  the  wave  energy. 

A  few  options  exist  for  implementing  diffusion  in  SWAN-X: 

•  use  of  a  weighted  mix  of  BSBT  and  SLl  sehemes, 

•  use  of  an  explieit  diffusion  seheme  that  ealeulates  diffusion  from  Eq.  (64)  using  information  at  the 
previous  time  level,  or 

•  use  of  an  implieit  diffusion  seheme  that  uses  Eq.  (64)  and  information  at  the  eurrent  time  level. 

6.1  Use  of  Numerical  Diffusion  (BSBT) 

This  option  would  involve  solving  the  propagation  equations  using  both  the  BSBT  seheme  and  the  SLl 
seheme.  The  averaging  of  the  two  sehemes  would  be  governed  by  a  weighting  faetor  that  would  be  ehosen 
based  on  eonsideration  of  the  desired  diffusion  level  (D.D,  and  Z)  )  as  well  as  the  relative  levels  of 
numerieal  diffusion  produeed  by  the  two  sehemes.  This  would  be  fairly  simple  to  implement,  but  the  diffu¬ 
sion  of  Eq.  (64)  would  be  impossible  to  reproduee  aeeurately  using  numerieal  diffusion.  Indeed,  the  repre¬ 
sentation  eould  be  quite  poor. 

6.2  An  Explicit  Diffusion  Scheme 

A  finite  differenee  representation  of  Eq.  (64)  eould  be  added  to  the  SLl  finite  differeneing  in  the  model. 
If  this  finite  differenee  seheme  for  diffusion  is  explieit  (based  on  independent  variables  ealeulated  at  the 
previous  time  step),  then  the  modifieation  to  the  model  would  be  minor.  This  alternative  has  the  disadvan¬ 
tages  with  regard  to  aeeuraey  and  stability:  fiilly  explieit  sehemes  tend  to  be  inaeeurate  eompared  to  sehemes 
that  use  information  at  the  eurrent  time  level  (e.g.,  eentered  fiilly  implieit  sehemes);  also,  fiilly  explieit 
diffusion  sehemes  are  invariably  eonditionally  stable.  Several  explieit  diffusion  sehemes  were  tested  outside 
the  model.  One  sueh  seheme  is  given  here.  It  was  derived  using  Taylor  series  with  the  objeetive  of  using  few 
points  that  are  outside  the  numerieal  steneil  used  by  the  SLl  seheme  (for  the  sake  of  eonvenienee): 


A 


(66) 


D 


-20AZ”-; +11A^-;_.) 


2D 


Here,  the  diffusion  eoeffieients  for  eaeh  wave  eomponent  are  assumed  uniform  in  spaee;  this  assump¬ 
tion  would  not  neeessarily  be  made  in  the  aetual  model.  Equation  (66)  was  (empirieally)  found  to  be  some¬ 
what  more  aeeurate  and  stable  than  a  setup  that  uses  traditional  lower  order  formulations  for  the  ^^Ydx^ 

d^N/ 

/dy^ 


terms. 
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6.3  An  Implicit  Diffusion  Scheme 

The  third  alternative  would  be  the  addition  of  an  implieit  seheme.  The  four-quadrant  algorithm  used  by 
SWAN  requires  that  all  eomputational  points  used  at  the  eurrent  time  level  must  be  upwind.  Unfortunately, 
no  uneonditionally  stable,  upwind,  implieit  diffusion  sehemes  are  known.  Therefore,  this  option  must  be 
ruled  out. 

7.  CONCLUSIONS 

The  higher  order  seheme  “SLl,”  taken  from  Stelling  and  Leendertse  (1992),  was  sueeessfiilly  imple¬ 
mented  in  an  experimental  version  of  SWAN  (“SWAN-X”).  The  modified  model  has  demonstrated  greatly 
redueed  diffusion  in  eases  where  diffusion  presents  problems  for  the  original  model.  However,  for  many 
nearshore  simulations,  SWAN  with  the  original  (BSBT)  seheme  is  more  appropriate  than  is  SWAN-X  with 
the  higher-order  SLl  seheme. 

Three  primary  faetors  support  this  argument: 

1)  The  SLl -based  model  is  relatively  slow  in  stationary  mode,  beeause  the  model’s  system  of  equa¬ 
tions  eannot  be  solved  by  simply  sweeping  through  the  geographie  grid  points.  Thus,  the  lower- 
order  seheme  should  be  used  in  stationary  simulations  for  whieh  the  higher-order  seheme  does  not 
offer  signifieant  benefit.  Most  small-seale  simulations  are  eondueted  in  stationary  mode. 

2)  The  original  SWAN,  whieh  uses  the  BSBT  seheme,  is  relatively  inexpensive  when  used  in  the 
nearshore  beeause  fine  geographie  resolutions  do  not  dietate  fine  temporal  resolutions.  SWAN-X 
with  the  SLl  seheme,  though  uneonditionally  stable,  does  not  have  this  freedom  from  eonsider- 
ations  of  Courant  number.  Thus,  when  fine  geographie  resolution  is  required,  smaller  time  step 
inerements  make  the  SLl -based  model  more  expensive  than  the  BSBT-based  model.  The  time  step 
size  must  be  eonsidered  for  the  SLl -based  model  in  both  stationary  and  nonstationary  mode. 

3)  The  diffusion  reduetion  of  the  SLl -based  model  is  insignifieant  for  many  simulations  on  the  eonti- 
nental  shelf,  beeause  numerieal  diffusion  at  these  seales  tends  to  be  small.  The  greatest  benefit  of 
the  higher-order  seheme  would  be  felt  more  at  oeeanie  seales.  There  are  some  notable  exeeptions  to 
this,  e.g.,  in  the  lee  of  islands,  where  diffusion  is  often  a  major  eoneern. 

Thus,  in  over-the-shelf  simulations,  the  higher-order  seheme  ean  be  used  when  numerieal  diffusion 
warrants  it.  This  greater  aeeuraey  eomes  with  a  signifieant  eost  when  the  simulation  is  either  stationary,  or 
high  resolution,  or  both.  Fortunately,  it  is  relatively  easy  to  prediet  whether  a  simulation  will  exhibit  high 
levels  of  diffusion,  so  a  user  ean  deeide  whether  the  potentially  high  eomputational  eost  is  warranted. 

Numerieal  diffusion  of  a  geographie  propagation  seheme  is  eaused  by  gradients  of  wave  aetion  aeross 
geographie  grid  points.  The  following  faetors  lead  to  sueh  gradients: 

•  nonstationary  input, 

•  nonuniform  input, 

•  obstaeles  in  the  model  domain,  and 

•  refraetion-indueed  features  (whieh  are  more  severe  at  loeations  with  rugged  bathymetry  or  spatially 
varying  eurrents). 

Also  note  that  broad  wave  speetra  will  tend  make  diffusion  less  notieeable  in  a  waveheight  field. 

At  oeeanie  seales,  simulations  are  nonstationary  and  relatively  nonuniform,  making  numerieal  diffu¬ 
sion  a  real  eoneern.  Nonstationary,  low-resolution  (Ax  >10  km)  simulations  ean  be  run  with  the  higher-order 
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scheme  with  very  little  added  computational  cost.  Thus,  the  implementation  of  the  higher-order  scheme  as 
an  option  in  the  official  version  of  the  model  is  one  major  step  toward  making  SWAN  a  viable  ocean-scale 
wave  model. 
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Appendix  A 

SWAN-X:  CHANGES  MADE  TO  THE  MODEL 


Al.  CHANGES  TO  CODE 

The  following  changes  were  made  to  the  code,  organized  by  file  and  subroutine.  Irrelevant  changes 
such  as  comment  statements  are  not  listed. 

swanprel.for: 

SUBROUTINE  SWREAD 

•  The  code  is  altered  such  that  the  variable  “ACl”  is  always  used  in  “STEL”  mode  (i.e.,  when  SLl 
scheme  is  used). 

•  The  “SCHEME”  command  is  added  to  user  input  (see  A2. 1  below). 

•  The  “COMPUTE”  command  is  modified  to  read  only  At  when  in  stationary  mode  with  SL 1  scheme. 

swanmain.for: 

SUBROUTINE  SWMAIN 

•  There  are  9  points  in  the  numerical  stencil  when  the  SLl  scheme  is  used. 

•  The  code  is  altered  such  that  the  variable  “AC2”  is  always  rolled  back  to  “ACl”  when  in  “STEL” 
mode. 

swancoml.for: 

SUBROUTINE  SWCOMP 

•  In  stationary  SLl  mode,  the  wave  action  variable  is  rolled  back  after  every  iteration  by  calling  the 
new  subroutine  ST  SNEXTI 

SUBROUTINE  SWOMPU 

•  Group  velocity  at  previous  time  level  is  added  to  the  variable  list. 

•  Necessary  points  are  added  to  numerical  stencil  if  SLl  scheme  will  be  used  at  given  grid  location. 

•  If  in  nonstationary  SLl  mode,  calculate  group  velocity  at  previous  time  step  (unlike  action  variable, 
group  velocity  is  not  saved  in  memory,  so  must  be  calculated). 

•  If  in  stationary  SLl  mode,  set  group  velocity  at  previous  time  step  equal  to  group  velocity  at  current 
time  step. 

•  If  in  SLl  mode,  subroutine  SPREDT  is  not  called.  This  modification  is  not  strictly  necessary. 
SUBROUTINE  ACTION 

•  If  in  SLl  mode  and  not  near  grid  boundary,  call  SANDL  instead  of  STRSXY. 

swancomS.for: 

SUBROUTINE  STRSXY 

•  When  in  “STEL”  mode,  action  variable  at  explicit  time  level  is  always  “ACl.” 

SUBROUTINE  SANDL  (new) 

•  This  new  subroutine  performs  the  same  function  as  subroutine  STRSXY,  but  calculates  gradients 
using  SLl  scheme. 
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A2.  CHANGES  TO  USER  INPUT 

A  small  number  of  changes  were  made  to  user  input  (provided  to  the  model  by  the  user  in  the  file 
“INPUT”): 

A2.1  “SCHEME”  Command 


I— >BSBT 

SCHEME  <  > 

I  STEL 


With  this  option,  the  user  specifies  the  numerical  scheme  to  be  used  for  geographic  propagation.  “BSBT” 
indicates  that  the  classic  first-order  upwind  implicit  scheme  will  be  used.  “STEL”  indicates  that  a  higher- 
order,  cyclic  scheme  will  be  used — the  =  0,  q  =  1/6”  scheme  of  Stelling  and  Leendertse  (1992).  (Note: 
the  BSBT  scheme  will  be  used  at  grid  boundaries,  regardless  of  this  setting.)  The  STEL  scheme  is  much  less 
diflusive  than  the  BSBT  scheme.  Thus,  it  should  be  used  in  cases  where  the  wave  field  is  highly  non- 
uniform  and/or  when  the  model  input  is  highly  nonstationary.  The  STEL  scheme  is  unconditionally  stable 
and  has  been  used  successfully  with  Courant  numbers  higher  than  50.  However,  the  scheme  often  produces 
oscillatory  behavior  at  high  Courant  numbers  (CFL»1),  so  it  should  be  used  with  caution. 

A2.2  “COMPUTE”  Command 


I— >SEC  I 

COMPUTE  ([deltc]  <  MIN  >) 

I  hr  I 

I  DA  I 


This  change  applies  only  to  the  model  that  uses  “MODE  STAT”  and  “SCHEME  STEL.”  With  this 
command,  a  time  step  size  is  specified.  The  begin  and  end  times  are  unnecessary;  the  number  of  time  step 
increments  is  specified  as  “itermax”  in  the  command  “NUM  ACCUR.” 

A2.3  “NUM  ACCUR”  Command 


NUMERIC  ACCUR  Idrell  [dabs]  [dtabs]  [npnts]  [itermax] 


Specification  of  this  command  is  unchanged.  However,  it  should  be  noted  that  when  in  stationary  mode 
while  using  the  “STEL”  scheme,  the  variable  “itermax”  refers  to  the  number  of  time  steps  that  the  model 
will  allow.  Time  steps  are  treated  as  iterations:  at  every  time  step,  convergence  is  checked  to  determine 
whether  the  model  has  reached  steady  state.  If  the  user  desires  multiple  iterations  at  each  time  step, 
nonstationary  mode  should  be  used  (note,  however,  that  convergence  is  not  checked  in  nonstationary  mode). 

A2.4  Caveats 

When  using  “MODE  STAT”  and  “SCHEME  STEL,”  the  two  caveats  described  below  apply. 
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Improving  the  Numerics  of  a  Third-Generation  Wave  Action  Model 


69 


A2.4.1  Caveat#! 

In  cases  where  very  small  time  increments  are  used,  the  model  solution  may  change  very  little  between 
time  steps;  there  is  a  danger  that  the  convergence  criterion  will  pass  as  a  result,  causing  the  model  to  end 
prematurely.  To  avoid  this,  the  user  may  specify  “drel,”  “dabs,”  and/or  “dtabs”  to  small  values,  or  specify 
“npnts”  to  a  high  value  (or  set  npnts  >  100  to  make  convergence  impossible). 

A2.4.2  Caveat  #2 

When  choosing  “itermax”  and  the  time  step  size  At  (“deltc”),  the  user  should  be  cognizant  of  the  result¬ 
ing  model  duration  —  the  wave  energy  should  be  given  time  to  fiilly  propagate  and  wind  generation  given 
time  to  produce  a  fully  arisen  sea  condition. 


Appendix  B 

DIFFRACTION  IN  SWAN 


Bl.  INTRODUCTION 

As  mentioned  in  Seetion  2.2,  the  possibility  of  modifying  SWAN  to  inelude  diffraetion  explieitly  is 
eurrently  being  eonsidered.  The  questions  eentral  to  the  debate  would  be  “when  does  diffraetion  make  a 
large  differenee?”  and  “are  these  situations  eommon  enough  to  warrant  modifieation  of  SWAN?”  The  most 
obvious  examples  of  diffraetion  are  assoeiated  with  struetures  sueh  as  jetties,  groins,  and  breakwaters.  In 
these  eases,  it  is  elear  that  diffraetion  will  move  energy  into  the  “shadow”  region  behind  the  strueture.  As 
there  are  no  other  physieal  proeesses  that  do  this,  elearly  diffraetion  is  important  here.  However,  one  might 
argue  that  simulations  with  struetures  sueh  as  this  tend  to  be  at  small  seale,  where  more  rigorous  models 
sueh  as  REF/DIF  (Kirby  and  Ozkan  1994)  and  CGWAVE  (Panehang  et  al.  1991) — ^both  of  whieh  inelude 
diffraetion  implieitly — ean  be  used  instead  of  SWAN. 

More  eommon  seenarios,  sueh  as  a  ease  of  wave  transformation  on  a  nonplanar  beaeh,  are  also  affeeted 
by  diffraetion,  as  diffraetion  tends  to  smooth  alongwave  variations  in  waveheight,  and  ereate  distinetive 
features  in  the  waveheight  field  in  loeations  where  alongwave  variation  is  espeeially  severe.  However,  the 
impaet  of  diffusion  in  sueh  eases  is  sometimes  not  notieeable,  due  to  multidireetional  wave  eonditions 
“blurring”  the  effeets  of  diffraetion.  The  sensitivity  of  diffraetion  effeets  to  wave  eonditions  is  not  obvious. 
Therefore,  a  test  ease  was  used  to  determine  this  sensitivity.  A  lab  experiment  was  available  that  was  well 
suited  for  this:  the  “WET  shoal”  of  Vineent  and  Briggs  (1989). 

B2.  EXPERIMENT  DESCRIPTION 

This  Vineent  and  Briggs  experiment  is  a  test  ease  of  wave  transformation  over  an  elliptieal  mound  (with 
a  eonstant  1.5 -ft  depth  over  the  rest  of  the  domain).  A  variety  of  irregular  wave  eonditions  were  used.  Input 
wave  speetra  were  well  doeumented,  as  were  waveheights  aeross  the  interior  of  the  domain  (espeeially 
within  the  shoal’s  region  of  infiuenee).  Figure  Bl  shows  the  bathymetry  and  gauge  loeations.  The  “non¬ 
breaking”  series  of  the  experiment  were  used  for  eomparison  to  SWAN.  A  waveheight  of  2.54  em,  a  wave 
period  of  1 .3  s,  and  an  ineident  wave  direetion  of  0°  were  used  for  all  eases  in  this  series;  Table  B 1  shows  the 
speetral  wave  eharaeteristies. 

In  this  experiment,  refraetion  ereates  a  region  of  energy  foeusing  direetly  behind  the  shoal.  Diffraetion 
aets  on  this  refraetion-indueed  feature.  To  demonstrate  the  eapability  of  a  phase-resolving  model  to  simulate 
diffraetion  implieitly,  the  monoehromatie  ease  (M2)  and  near-monoehromatie  ease  (U4)  are  eompared  to 
REF/DIFl  (Kirby  and  Dalrymple  1993).  The  model  does  an  extraordinary  job  of  eapturing  the  two  diffrae- 
tion-indueed  features  adjaeent  to  the  eentral  refraetion-indueed  feature.  The  absenee  of  diffraetion  is  appar¬ 
ent  in  the  SWAN  result  for  the  near-monoehromatie  ease  (Fig.  B2).  However,  the  model  seems  to  suffer 
from  the  omission  of  diffraetion  only  for  the  least  realistie,  monoehromatie  eases  (Fig.  B3). 
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Fig.  B1  —  Bathymetry  and  gauge  loeations  for  the  Vineent  and  Briggs  shoal  experiment 


Table  B1  —  Spectral  Characteristics  of  the  Nonbreaking  Series  of 
Vincent  and  Briggs  (1989) 


Case  ID 

Directional  Spectrum 

Frequency  Spectrum 

M2 

Unidirectional 

Monochromatic 

U3 

Unidirectional 

Broad 

N3 

Narrow 

Broad 

B3 

Broad 

Broad 

U4 

Unidirectional 

Narrow 

N4 

Narrow 

Narrow 

B4 

Broad 

Narrow 
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Fig.  B2  —  Comparison  of  normalized  waveheights  along  a  lateral  transeet  behind 
the  shoal.  REF/DIFl  output  is  eompared  to  the  monoehromatie  ease  and  SWAN 
output  to  the  near-monoehromatie  ease. 


distance  along  wavemaker  (m) 

Fig.  B3  —  Comparison  of  normalized  waveheights  along  a  lateral  transeet 
behind  the  shoal.  SWAN  output  for  the  six  of  the  nonbreaking  wave  eondi- 
tions  are  eompared  to  eorresponding  data. 
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B3.  CONCLUSION 

This  is  merely  a  single,  simple  experiment  of  diffraetion  on  a  nonplanar  bathymetry.  However,  it  sug¬ 
gests  that  under  realistie,  multidireetional  wave  eonditions,  and  in  the  absenee  of  wave-bloeking  struetures, 
diffraetion  has  only  a  minor  impaet  on  gross  energy  estimates. 

When  information  is  needed  beyond  simple  energy  estimates,  the  inelusion  of  diffraetion  may  be  a 
serious  eonsideration.  One  would  expeet  the  impaet  of  diffraetion  on  direetional  distributions  of  wave  en¬ 
ergy  to  be  nontrivial  in  many  eases.  Also,  if  phase-related  information  is  needed,  e.g.,  wave  asymmetry,  the 
use  of  a  phase-resolving  model  that  ineludes  diffraetion  is  elearly  neeessary. 
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Appendix  C 

RECENT  DEVELOPMENTS  —  STATIONARY  HIGHER-ORDER  SCHEMES 


Cl.  INTRODUCTION 

Although  a  pseudo-stationary  mode  for  SWAN-X  w/SLl  was  sueeessfully  implemented  (Seetion  4.3.1), 
it  is  elear  that  the  SLl  seheme  is  not  suited  for  stationary  eomputations.  Due  to  the  downwind  point  in  its 
numerieal  steneil,  the  seheme,  operating  at  a  single  time  level,  eannot  simply  sweep  through  the  geographie 
grid.  The  pseudo-stationary  method  of  time-stepping  until  steady  state  is  impraetieal  due  to  the  large  number 
of  time  steps  typieally  required.  Iterative  methods  for  a  stationary  SLl  seheme  have  proven  unreliable. 
Thus,  we  have  endeavored  to  find  a  seheme  more  suited  for  stationary  ealeulations. 

C2.  TESTS  OF  STATIONARY  SCHEMES  OUTSIDE  OF  SWAN 

Beeause  we  wanted  to  avoid  the  use  of  iterative  methods,  we  imposed  a  requirement  on  the  new  seheme: 
it  would  need  to  be  entirely  upwind.  We  eompared  the  results  of  simple  models  based  on  various  upwind 
sehemes  to  exaet  solutions  for  the  ease  of  a  submerged  breakwater  with  gap  (sinee  we  were  not  interested  in 
nonstationary  ealeulations,  a  large  set  of  test  eases  was  not  neeessary).  Tests  were  similar  in  setup  to  those 
shown  in  Figs.  1  and  24  of  the  body  of  this  report.  Results  were  eompared  to  the  exaet  solutions  along  the 
“high-x”  open  boundary. 

Some  of  the  more  notable  eontenders: 

1)  the  implieit  half  (i.e.,  the  upwind  stage)  of  the  SLl  seheme, 

2)  a  seheme  derived  from  Taylor  Series,  with  three  upwind  points  in  the  steneil, 

3)  the  Box  seheme, 

4)  the  Seeond  ORder  UPwind  Implieit  seheme  (Seetion  3.1 .2.2),  “SORDUPI,” 

5)  the  one-dimensional  NISL  seheme  (Seetion  3.1 .2.10),  made  two-dimensional  and  stationary  by 
replaeing  the  time  dimension  with  a  seeond  spaee  dimension,  and 

6)  same  as  (5),  but  using  the  SLl  seheme. 

Of  sehemes  (1)  through  (4),  seheme  (4)  displayed  the  least  severe  numerieal  oseillations  and  the  most 
severe  numerieal  diffusion  (see  Fig.  Cl).  As  numerieal  oseillations  were  deemed  to  be  more  of  a  eoneern 
than  diffusion  (given  the  faet  that  all  four  were  mueh  less  diffusive  than  the  BSBT  seheme),  seheme  (4)  was 
judged  the  most  favorable  of  these  four  sehemes. 

Of  the  two  “traded”  sehemes  (5  and  6),  the  NISL  seheme  is  elearly  more  aeeurate  (e.g..  Fig.  C2);  in  faet, 
it  was  the  most  aeeurate  seheme  tested.  However,  there  are  potential  diffieulties  assoeiated  with  the  NISL 
seheme: 
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Fig.  Cl  —  Result  from  several  stationary  models  of  the  submerged  breakwater  with  gap 
test  ease,  eompared  to  exaet  solution.  Six  propagation  angles  were  tested;  the  result  for 
theta  =  35  degrees  is  shown  here.  Distanee  is  nondimensionalized  by  Courant  number. 
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Fig.  C2  —  Result  from  two  stationary  models  of  the  submerged  breakwater  with  gap 
test  ease,  eompared  to  exaet  solution.  These  two  models  are  based  on  the  two  “traded” 
numerieal  seheme.  Six  propagation  angles  were  tested;  the  result  for  theta  =  35  de¬ 
grees  is  shown  here.  Distanee  is  nondimensionalized  by  Courant  number. 
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a)  The  variable  numerieal  steneil  of  the  NISL  seheme  is  dependent  on  the  “traded”  Courant  number 
(e.g.,  CFL  =  (C^^Ay)/(C^^Ax)).  The  subroutine  strueture  of  SWAN  does  not  presently  allow  for  this, 
so  nontrivial  ehanges  to  the  eode  would  be  neeessary.  Also,  at  a  given  geographie  grid  loeation,  the 
numerieal  steneil  for  geographie  propagation  varies  in  speetral  spaee.  NISL  would  require  eonsider- 
ably  more  RAM,  due  to  the  need  to  store  and  for  every  steneil  (note  that  with  the  dynamie 
memory  of  Fortran90,  this  would  be  less  of  a  eoneem). 

b)  It  is  uneertain  how  well  NISL  would  perform  for  large  Courant  numbers  when  varies  along  the 
eharaeteristie  used  by  this  semi-Lagrangian  seheme.  Numerieal  experiments  suggest  serious  prob¬ 
lems  with  mass  eonservation  for  sueh  eases. 

In  eontrast,  the  SORDUPI  seheme  appeared  simple  to  implement  in  SWAN,  with  no  potential  draw- 
baeks.  Therefore,  a  seeond  experimental  version  of  SWAN  (“SWAN-XX”)  was  ereated  that  uses  the  SORDUPI 
seheme.  Implementation  of  NISL  in  SWAN  is  still  a  possibility,  but  has  not  yet  been  undertaken. 

As  of  this  writing,  there  are  four  sehemes  that  are  either  already  in  SWAN  or  still  being  eonsidered  for 
inelusion  in  SWAN.  Based  on  the  breakwater  tests  deseribed  above,  the  four  sehemes  were  seored  based  on 
average  amplitude  preservation  (note  that  this  seoring  does  not  take  phase  error  into  aeeount)  as  shown  in 
Tabled. 


Table  Cl  —  Stationary  Breakwater  with  Gap  Test:  Amplitude  Preservation 


Scheme 

Amplitude  Preservation  (%) 

BSBT 

28 

SORDUPI 

55 

SLI  (nonstationary;  with  a  small  time  step) 

69 

NISL 

76 

C3.  SOUTHERN  CALIFORNIA  BIGHT 

The  SORDUPI-based  model  was  applied  to  the  Southern  California  Bight  test  ease  (Seetion  5.8).  Com¬ 
putation  times  for  this  model  were  eneouraging.  The  SORDUPI  seheme  uses  the  same  simple  sweeping 
teehnique  as  that  used  by  the  BSBT-based  model.  No  iterations  are  required  for  geographie  propagation 
(although  some  are  required  for  speetral  propagation).  The  additional  math  operations  used  by  the  seeond- 
order  seheme  inereased  eomputation  time  by  approximately  10%  per  iteration.  However,  in  the  100-m- 
resolution  ease,  this  was  offset  by  the  slightly  lower  average  number  of  iterations  required  by  the  SORDUPI 
model.  Table  C2  eompares  the  total  eomputation  times  on  an  UltrabO  (360  MHz)  Sun  workstation  with  the 
different  sehemes  and  geographie  resolutions. 

Outputs  of  the  models  were  eompared  at  the  loeations  of  the  four  pressure  gauges.  Error  in  the  SWAN 
model  results  at  gauge  10  (the  gauge  in  the  lee  of  San  Miguel  Island)  is  dominated  by  diffusion  effeets  and, 
therefore,  provides  a  good  indieation  of  relative  diffusivity  of  the  models.  The  results  are  as  expeeted:  at 
equivalent  resolution,  the  SLl  seheme  is  somewhat  less  diffusive  than  the  SORDUPI  seheme,  whieh  is  in 
turn  mueh  less  diffusive  than  the  BSBT  seheme  (Fig.  C3). 
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Table  C2  —  Southern  California  Bight:  SWAN  Computation  Times 


Scheme 

Resolution  (Ax  and  Aj) 

Total  computation  time^ 

BSBT 

100  m 

27  hours,  34  minutes 

SORDUPI 

100  m 

27  hours,  27  minutes 

SLl 

100  m 

Approximately  5  months  (this  set 
was  not  completed) 

BSBT 

200  m 

7  hours,  7  minutes 

SORDUPI 

200  m 

8  hours,  7  minutes 

SLl 

200  m 

Approximately  19  days 

BSBT 

400  m 

1  hour,  51  minutes 

SLl 

400  m 

55  hours,  15  minutes 

*  Durations  given  are  totals:  27  simulations  were  performed  with  eaeh  model 


Energy  Estimates  at  San  Miguel  Gage  (Gage  #10) 


Fig.  C3  —  Time  series  eomparisons  of  models  to  data  for  the  Southern  California  Bight  ease.  Here, 
SWAN  models  of  equivalent  geographie  resolutions  are  eompared. 
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When  models  of  approximately  equivalent  eomputation  time  are  eompared,  the  SORDUPI  seheme  is 
most  aeeurate  (Fig.  C4).  The  SLl-based  model  had  to  be  run  at  a  grid  resolution  of  Ax  =  Ay  =  400  m  in  order 
to  finish  within  a  reasonable  time  period.  At  this  resolution,  the  SLl-based  model  suffered  from  signifieant 
diffusion  (and  perhaps  other  resolution-related  effeets). 


Energy  Estimates  at  San  Miguel  Gage  (Gage  #10) 


Fig.  C4  —  Time  series  eomparisons  of  models  to  data  for  the  Southern  California  Bight  ease.  Here, 
SWAN  models  of  approximately  equivalent  eomputational  requirements  are  eompared. 


C4.  CONCLUSION 

The  eurrent  version  of  SWAN-X  has  two  sehemes  available  for  geographie  propagation:  the  BSBT 
seheme  (whieh  is  very  diffusive)  and  the  SLl  seheme  (whieh  is  not  suited  for  stationary  eomputations). 
Thus,  for  stationary,  diffiision-prone  simulations,  a  seeond  alternative  is  needed.  The  seeond-order  upwind 
implieit  (SORDUPI)  seheme  was  tested  and  results  look  promising. 

A  model  using  the  SORDUPI  seheme  would  be  mueh  less  diffusive  than  the  BSBT-based  model.  It 
would  be  somewhat  more  diffusive  than  the  SLl  model,  but  unlike  the  SLl  model,  is  eeonomieal.  In  faet,  it 
is  probably  no  more  expensive  than  the  BSBT  model.  Nonphysieal  oseillation  with  the  seheme  appears  to  be 
well  within  aeeeptable  levels.  Thus,  the  SORDUPI  model  is  a  mueh  better  ehoiee  for  a  large  majority  of 
stationary  eomputations  than  either  the  BSBT  seheme  or  the  SLl  seheme. 


